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1  Introduction 


We  give  below  a  rather  complete  account  of  the  integral-equation  formulations  as  they  are 
being  implemented  in  our  solver  for  elasticity  problems. 

We  start,  in  Sec.  3,  with  a  discussion  of  several  equivalent  forms  of  the  Lame  differen¬ 
tial  equation  in  elasticity,  including  also  its  first-order  representation  as  a  set  of  coupled 
equations  for  the  displacement  and  the  stress  tensor.  In  addition,  we  derive  forms  of  the 
differential  equations  with  separated  terms  describing  solution  in  the  “background  medium” 
(such  as  air)  and  interaction  terms  describing  effects  of  the  deviations  of  the  medium  prop¬ 
erties  from  those  of  the  background  material. 

Next,  we  discuss  (in  Sec.  4)  the  Green  functions  associated  with  the  Lame  equation  for 
the  displacement  field.  In  particular,  we  obtain  there  a  form  of  the  Green  function  which 
explicitly  exhibits  a  nonsingular  behavior  of  its  dyadic-derivative  part,  which  facilitates 
discretization  of  the  resulting  integral  equations. 

In  the  following  Sections,  5  and  6,  we  are  concerned  with  integral-equation  formulations 
in  elasticity,  in  three  cases:  (A)  purely  surface  (boundary-element)  equations,  (B)  purely 
volumetric  (Lippmann-Schwinger,  or  L-S)  equations,  and  (C)  a  coupled  system  of  volume 
and  surface  equations.  The  last  of  these  is  a  novel  approach  we  developed,  mostly  in  order 
to  be  able  to  efficiently  model  geometry  components  -  in  our  case,  the  middle  and  inner 
ear  -  characterized  by  small  sizes  and  intricate  surfaces,  and  embedded  in  a  larger  volume 
of  an  inhomogenous  material. 

2  Summary  of  the  integral  equations  in  elasticity 

We  briefly  the  describe  here  the  three  considered  types  of  the  integral  equations,  and  then 
list  the  genral  forms  of  the  equations  themselves. 

(A)  Surface  integral  equations.  The  surface  integral  equations  -  or  boundary  integral 
equations,  BIEs  -  are  applicable  to  piecewise  homogeneous  materials,  and  provide  solutions 
for  the  displacement  and  traction  fields  defined  on  interfaces  separating  different  material 
regions.  Fields  in  the  individual  regions  are  described  in  terms  of  the  appropriate  Green 
functions  for  elastic  materials.  We  envisage  using  this  type  of  equations  in  several  situations, 
such  as 

1.  Modeling  of  man-made  objects  of  possibly  complex  geometrical  shapes,  but  consisting 
of  only  few  homogeneous  materials. 

2.  Solution  of  the  auxiliary  surface  problems  arising  in  solution  of  the  volumetric  integral 
equations  for  materials  characterized  by  large-contrast  discontinuities  in  the  material 
properties  (we  return  to  this  problem  below). 

3.  More  generally,  verification  and  checks  of  the  solutions  obtained  with  the  volumetric- 
equations  code  in  the  case  of  piecewise  homogeneous  materials. 
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(B)  Volumetric  integral  equations.  Volumetric  integral  equations,  on  the  other  hand, 
can  be  used  for  inhomogeneous  media,  with  (generally)  different  material  properties  as¬ 
signed  to  the  individual  tetrahedra  into  which  the  volume  has  been  discretized.  We  use 
here  Lippnrann-Schwinger  equations  with  the  Green  function  associated  with  the  infinite 
(unbounded)  background  medium  -  in  our  case,  air. 

We  consider  two  types  of  volumetric  integral  equations: 

(i)  Equations  derived  from  differential  equations  in  their  first-order  form.  In  this  case 
the  unknowns  in  the  integral  equations  are  the  displacement  and  stress  tensor  fields 
defined  in  the  considered  volume. 

(ii)  Equations  derived  from  differential  equations  in  their  second-order  form.  In  this  for¬ 
mulation  the  unknowns  are  only  the  components  of  the  displacement  field. 

In  both  cases  we  discretize  the  integral  equations  in  terms  of  piecewise  linear  basis 
functions.  Each  such  basis  function  is  associated  with  the  vertices  of  the  tetrahedral  mesh, 
and  supported  on  sets  of  tetrahedra  adjacent  to  the  considered  vertex. 

In  designing  our  integral-equation  formulation  we  pay  a  particular  attention  to  the 
problem  of  a  possible  discontinuous  behavior  of  the  material  properties,  which  is,  clearly, 
always  present  when  considering  a  mechanically  dense  (biological)  material  immersed  in  air. 

Such  problems  are  known  to  cause  difficulties  in  solving  integral  equations  in  acoustics 
and,  in  that  case,  we  have  devised  an  approach  [1]  in  which  the  original  system  of  equations 
is  reformulated  in  terms  of  a  surface  problem  associated  with  the  contrast  interface(s) 
(characterized  by  a  large  ratio  of  densities  of  the  adjacent  materials)  and  a  “residual” 
volumetric  problem.  The  analogous  problem  in  elasticity  appears  to  be  considerably  more 
complex,  and  we  have  put  much  effort  into  deriving  appropriate  forms  of  integral  equations, 
in  both  first-  and  second-order  formulations  ((i)  and  (ii)  above). 

(C)  Coupled  volume-surface  integral  equations.  A  coupled  system  of  volumetric  and 
surface  integral  equations  arises,  e.g.,  when  a  homogeneous  material  region  is  embedded,  as 
an  inclusion,  in  an  inhomogeneous  material.  We  briefly  summarize  the  structure  of  such  a 
system  of  equations  below. 

General  forms  of  the  integral  equations.  We  give  below  formulae  for  the  systems  of 
integral  equations  used  in  this  Report.  The  derivations  are  given  and  properties  of  these 
equations  are  discussed  in  detail  in  Sections  5  and  6. 

(A)  Surface  integral  equations.  We  first  present  the  general  form  of  the  surface  in¬ 
tegral  equations  for  a  set  of  homogenous  regions  Qm  separated  by  interfaces;  one  of  these 
regions,  fig,  is  the  unbounded  background  medium.  The  displacement  and  traction  fields 
are  assumed  to  be  continuous  across  the  interfaces. 

The  resulting  system  of  integral  equations  simply  consists  of  two  equations  per  interface 
(oriented  surface)  Smn  separating  the  regions  Q,m  (on  the  negative  side  of  the  interface)  and 
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iln  (on  its  positive  side) 


5  u(r)  +  [  d2 r'  [rm(t,  r')  •  u(r')  +  (?£(  r,  r)  •  t(r')] 

Y  J  d  V  [r£(t,  r')  •  u(r')  +  G^(:  r,  r)  •  t(r')]  =  Sm0  um(r)  for  r  G  Smn  ,  (2.1a) 

s.  ean™ 

i^n 

\  u(r)  -  /  d2r'  [r^(t,  r')  •  u(r')  +  G^(r,  r')  •  t(r')] 

J  S rnn 

+  Y  j  d2r/[r^(t,r/)-u(r')  +  G^(r,r/)-t(r/)]  =5n0um(r)  for  r  €  Smn  .  (2.1b) 

S  GdClr,  °nj 

nj  n 

j^m 

With  reference  to  Fig.  1,  Eq.(2.1a)  represents  contributions  to  the  displacement  field  u  on 
the  interface  Smn  due  to  the  displacement  and  traction  fields  u  and  t  on  the  same  interface 
(the  first  integral)  and  on  other  interfaces,  Sim,  forming  boundaries  of  the  region  VLm  with 
other  regions  VLi:  i  ^  n.  These  intervals  involve  Green  functions  Gm  and  Fm  (defined  by 
Eqs.  (4.20)  and  (5.4)  in  Sec.  5.1),  describing  propagation  of  the  fields  in  the  region  f 1m. 
Similarly,  Eq.(2.1b)  represents  contributions  to  the  field  u  on  the  interface  Smn  due  to 
the  fields  on  the  boundaries  of  the  other  region,  Fln,  adjacent  to  the  interface.  The  r.h.s.s 
of  Eqs.  (2.1)  are  the  incident  fields  due  to  distant  sources  in  the  region  (hence  the 
delta-functions  Sm0  and  <5n0 ) . 


Figure  1:  A  schematic  representation  of  regions  f 1  and  interfaces  S  appearing  in  integral 
equations  2.1. 

(B)  Volumetric  integral  equations.  We  obtained  our  volumetric  equations  are  the 
Lippmann-Schwinger  (L-S)  equations,  describing  the  elastic  fields  as  propagating  in  the 
background  medium  (characterized  by  some  density  p0  and  Lame  coefficients  A0  and  //0), 
modified  by  the  presence  of  the  actual  medium  with  material  parameters  p.  A,  and  p.  The 
interaction  of  the  elastic  wave  with  the  medium  may  be  described  in  many  equivalent  ways, 
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resulting  in  various  forms  of  the  L-S  equation.  Our  forms  of  these  equations  are  not  the 
conventionally  used  ones,  as  we  took  special  care  to  ensure  their  favorable  properties  in 
problems  involving  large  density  contrasts. 


(B  i)  Equations  in  the  “first-order”  form.  In  this  formulation,  derived  from  the 
Lame  equation  in  its  first-order  form,  the  unknowns  are  the  velocity  Uj  :=  — i  kui  (where 
k  is  the  wave  number  in  the  background  medium),  the  pressure  p  :=  —  and  the 

symmetric  traceless  part  atJ  of  the  stress  tensor  Ttj.  The  equations  involve  the  Green 
function  g(r)  =  exp(i/cr) /(47rr)  of  the  Helmholtz  equation  in  the  background  medium,  as 
well  as  its  derivatives, 

5mn(r)  :=  (|  5mn  k2  +  dmdn)  g{ r)  .  (2.2) 

The  obtained  system  of  integral  equations,  Eqs.  (6.35),  is  then 

—  *>i(r)  +  [  d  V  (dtdm  g(r  -  r'))  ~  l)  vm{r') 

Po  J  VP  o  / 

+  ^2  3m  (di  VW  +  dm  ^(r)) 

+  ^2  j  dV  ( 9i  9mn (r  -  r0)  d'm.  «n(r') 

+  YoJ  d3r>  ( di  p(r  -  r/))  (<p(1-/)  _  x)  p(r/) 

=  ^n(r)  ,  (2.3) 

p{ r)  —  k2  J  d3/  g( r  —  r')  (</?(r')  —  l)  p( r') 

+  ^  ^  dm  Vm(r)  +  \  j  dV  Pmn(r  “  r0  d'm  Gx(r0 

=  pin(r)  , 

°ij(r)  -  J,  [ 9i  vj(r)  +  dj  vi(r)  -  |  6ij  dm  vm(r)] 

=  O',™  (r)  , 

with  a  dimensionless  material  parameter 


V  = 


A  +  §/j 


(2.4) 
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(B  ii)  Equation  in  the  “second-order”  form.  In  this  case  we  have  only  a  single 
equation  Eq.(6.51)  for  the  three-component  displacement  field  u, 


«i(r)  -  k  2  dj  (1  -  £A(r))  Si:jdkuk( r)  -  ^(r)  (^(r)  +  djUi( 


-  k 


-2 


8,; 


Po 

P(r) 


[V\(r)  5ij  dkuk(r)  +  V„(r)  [diUji r)  +  djUi{ 

—  k~2  J  d 3r'  di  dj  g( r  —  r') 

(!  -  £A(r'))  Stj  d'kuk(r')  -  ^(r')  (d[uj{r')  +  8ju,( r' 

-  fc"2  J  dV  9-  <?(r  -  r')  ^ 

7?A(r')  ^^(rO  +  ^r')  (3,/ui(r/)  +  9'u,(r/ 

=  , 

with  the  auxiliary  material-dependent  coefficients  defined  as 

t  =  PoA  r,  =  A  r,  =  ii 

^>A  „  \  ’  ‘’/i  n  \  ’  1 A  \  ’  m  \  • 

P  Ao  P  Ao  Ao  Ao 


(2.5) 


(2.6) 


Properties  of  the  integral  equations  in  high-contrast  problems.  The  volumetric 
integral  equations  (2.3)  and  (2.5)  have  quite  different  (although  equivalent)  forms;  they 
share,  however,  common  features,  which  become  relevant  in  problems  involving  large  con¬ 
trast.  Such  cases,  in  our  applications,  are  characterized  by  large  ratios  of  material  density 
and  the  Lame  parameter  A  values  in  the  considered  material  and  in  the  background  medium, 
with  moderate  values  of  the  wave  propagation  speed,  i.e., 


P_ 

Po 


A 


~  —  3>  1  . 
An 


(2.7) 


In  this  limit  only  some  terms  in  the  integral  equations  are  dominant  and,  moreover,  they 
represent  contributions  of  interfaces  at  which  there  occur  large  jumps  in  the  parameters 
p  and  A.  This  structure  of  the  equations  facilitates  their  discretization  and  solution  in 
high-contrast  problems. 


(C)  Coupled  volumetric  and  surface  integral  equations.  A  simple  example  of  a 
system  involving  both  volumetric  and  surface  fields  is  visualized  in  Fig.  2.  In  this  case 
surface  fields  u  and  t  are  defined  on  the  boundary  dVtm  of  a  homogeneous  region  Qm 
embedded  in  an  inhomogeneous  region  fi,  supporting  the  volumetric  field  u. 
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Figure  2:  A  schematic  representation  of  coupled  integral  equations:  a  volume  integral  equa¬ 
tion  for  the  displacement  field  u  in  the  inhomogeneous  region  12  and  a  surface  integral 
equation  for  the  displacement  and  traction  fields,  u  and  t,  on  the  boundary  of  the  homoge¬ 
neous  material  region  12m  embedded  in  12. 

In  notation  similar  to  that  in  Eqs.  (2.1),  the  resulting  coupled  integral  equations  have 
then  the  general  form 


for  r  €  12  , 


for  r  e  <9!2m  , 


for  r  G  <912m  , 


(2.8a) 

(2.8b) 

(2.8c) 


where  <S(u(r))  is  a  functional  of  the  displacement  field,  representing  a  volumetric  source, 
whose  specific  form  depends  on  the  implementation  of  the  volumetric  (Lippmann-Schwinger) 
equations.  The  Green  functions  labeled  with  the  index  m  refer  to  the  material  of  the  region 
12m,  while  those  without  the  label  are  the  background-medium  Green  functions. 
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3  Differential  equations  in  elasticity 


3.1  Second-order  differential  equations  for  isotropic  media 

General  Lame  equation.  From  the  frequency-domain  elasticity  equation  for  the  dis¬ 
placement  field  u(r), 

u2  pUt  +  djT^  =  0  , 


the  constitutive  relation 


Tij  —  Cijkl  ekl  ’ 

and  the  definition  of  the  strain  tensor, 

eij  =  \  (di  uj  +  dj  ui )  > 

one  obtains,  upon  setting 

Cijki  =  ^  $ij  3 kl  +  P  {&ik  3jl  +  ^ a  3jk)  t 
the  basic  form  of  the  Lame  equation 

io2  pui  +  di{\  dj  Uj)  +  dj  [n  [dt  Uj  +  dj  u^j]  =0 


(3-1) 

(3.2) 

(3-3) 

(3-4) 

(3-5) 


The  above  equation  is  equivalent  to  the  set  of  two  first  order  equations,  Eq.(3.1)  and  the 
definition  of  the  stress  tensor  in  terms  of  the  displacement, 


Tij  =  A  sij  dk  uk  +  M  {di  Uj  +  dj  ut)  .  (3.6) 

In  order  to  derive  the  L-S  equation  for  elasticity,  we  now  represent  Eq.(3.5)  in  a  form 
exhibiting  the  background-medium  operator  and  interaction  terms.  The  background  mate¬ 
rial  is  characterized  by  the  density  p0  and  the  Lame  coefficients  A0  and  /i0  =  0,  such  that 
the  wave  vector  k  in  that  medium  is 


k2  =  ^  cu2 


In  terms  of  k.  the  basic  form  of  the  Lame  equation  (3.5)  is 

f  (a,  »j  +  @j «.) 


k2  —  ^  +  d.:  (  -^-  dj  Uj  |  +  dj 


Po 


L'H) 


=  0 


(3.7) 


(3.8) 
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Lame  equation  in  L-S  form.  We  will  use  three  alternative  forms  of  the  differential 
equation,  equivalent  to  Eq.(3.8), 


(d?dJ  +  Sij  k2)  Uj  -  k2  (*  -  £)  ui  -  di  (X  -  diu3 

+  93  Yn  (9iU3  +  93Ui) 

( di°j  +  Sij  k2)  U3  ~  9i 


1-^1^ 


+  —  dj 
p  3 


9id3  +  Sij  k  U3  -  d3 


-  d 


p\ 


(>  -  Q  5ij  dlUl  -  N-  [diuj  +  d. 


=  0  , 

£  (w  +  dY 


-  1 1'.—  |  £  d  u  -  o  , 
'  P  )  Ao  3  3 


'3Ui 


Po 

> 


Vx  SijdiUi  +  Vp  [diUj  +  djU, 


=  0 


where  the  dimensionless  parameters  £  ,  r/A,  and  r)  are  defined  as 


(3.9a) 


(3.9b) 


(3.9c) 


PO  £  _  Po  I1 

p  X0  1  M  P  A0 


(3.10) 


and 

P\  =  Y'  ^  =  Y  ’  (3-n) 

Ao  Ao 

In  air  we  have,  by  definition,  =  r]x  =  1  and  £,  =  r//t  =  0,  while  in  typical  biological  media 
p0//0  ^  1)  V\  1)  and  r//(  3>  1,  but  the  coefficients  and  ^  (which  are  proportional  to 
the  refraction  coefficients  squared)  remain  of  order  1. 

The  heuristics  behind  Eq.(3.9b)  involves  representing  the  pressure  in  terms  of  the  dis¬ 
placement,  p  =  —XdjUi,  “removing”  one  differentiation,  and  adding  the  shear  term  with  the 
Lame  coefficient  p.  An  important  feature  of  Eq.(3.9b)  is  the  last  term  with  the  manifestly 
appearing  gradient  of  p0/p. 

Similar  remarks  apply  to  Eq.(3.9c),  in  which,  additionally,  the  terms  proportional  to  Y 
and  Y  combine  to  form  expressions  proportional  to  the  stress  tensor  (see  Eq.(6.52)  below). 


3.2  Second-order  differential  equations  for  anisotropic  media 

In  the  following  we  will  also  consider  equations  for  anisotropic  media;  the  main  reason 
for  this  is  that  we  will  impose  kinematic  constraints  in  elastic  shell  theories  by  means  of 
anisotropic  material  properties.  In  addition,  we  may  need  to  model  anisotropic  shells,  such 
as  the  basilar  membrane. 

For  an  anisotropic  medium  the  elasticity  tensor  C  of  Eq.(3.2)  is  a  general  fourth  rank 
tensor  satisfying  the  symmetry  relations 

Cijkl  =  Cjiki  >  ^ijkl  =  C ijlk  >  C ijkl  =  C klij  ’  (3-12) 

and  the  condition  that  the  quadratic  form 


W:=C 


ijkl  ^ ij  ^ kl 
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(3.13) 


is  positive-definite.  The  symmetry  properties  (3.12)  imply  that  the  tensor  C  involves  21 
independent  material  parameters. 

With  a  general  elasticity  tensor  C,  the  frequency-domain  form  of  Eq.(3.1)  becomes,  in 
analogy  to  Eq.(3.8), 

fc2  ^  Ui  +  2^  dj  \-Ciikl  ^ 9k  Ul  +  °l  Uk )]  =  °  ’  (3-14) 

3.3  First-order  differential  equations  for  isotropic  media 

General  first-order  equations.  The  second-order  Lame  equation  (3.8)  can  be  then 


written  as  a  system  of  two  first-order  equations 

k2  —  ui  +  d  r  - ■  =  0  ,  (3.15a) 

Tij  -  A  Sij  dkuk~P  ( di  uj  +  dj  ui)  =  0  >  (3.15b) 

containing  no  derivatives  of  material  parameters. 

As  a  matter  of  notation,  we  define  the  “velocity”  v  as 

v  =  —  i  k  u  (3.16) 

and  rewrite  the  above  equations  in  the  form 

ik-vi  +  -^-diTii  =  0,  (3.17a) 

Po  Ao 

i  k  +  A  dk  vk  +  n  (dl  Vj  +  dj  v^j  =  0  .  (3.17b) 


Next,  we  try  to  represent  these  equations  in  such  a  form  that  the  material  parameters  do 
not  appear  as  factors  of  the  differential  operators.  At  the  same  time,  it  should  be  possible 
to  split  those  equations  into  terms  describing  fields  in  the  background  medium  (p  =  p0, 
A  =  A0,  p  =  0)  and  interaction  terms;  this  structure  is  required  in  order  to  derive  the  L-S 
equations. 


A  straightforward  decomposition  into  background  and  interaction  terms.  We 

represent  now  Eqs.  (3.17)  as 


1  k  dk  +  hi  (y  _ d*  ** +  (d<  vi +aivi)=0 


JCF 

=  (V  +  V)F=  0 

T'wv  I  tavo; 

uik  i  Uikl 

i  k6ik\  6ikdi 

'TACJV  ! 

uijk  ;  uijkl 

5ijdk\[k5ik  5ji 

(3.18a) 

(3.18b) 

(3.19) 

(3.20) 
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(3.21) 


vvu; 

vikl 

ik5lk(X/\0-l) 

0 

v-I 

vijkl 

$ ij  (-V “  !)  &k  +  /V ao  (f>ik  dj  +  fijk 

)  0 

and 


1 

K 

1 _ 

i 

c2 

_ i 

i 

! _ i 

i 

_ 1 

i 

e 

i _ ! 

aij/\ 

(3.22) 


The  Green  function  corresponding  to  the  background-medium  differential  operator  (3.20) 
is  defined  by  the  equation 


vg{ r)  =  -X<J2(r)  , 


(3.23) 


where  X  is  the  unit  tensor  (we  discuss  more  general  Green  functions  for  elastic  media  in 
Sec.  4). 

By  taking  the  Fourier  transform  of  Eq.(3.23)  we  can  obtain  an  algebraic  equation  for 
the  Fourier  representation 

6(r)=/|^elqre(q)  (3.24) 

of  the  Green  function.  It  has  the  form 


'lkSik 

-i  sik<h 

~W 

y  km 

~VCd 

y  kmn 

~'l6ijQk 

[k6ik6ji 

a*  < 

~UJU) 

klmn 

(q) 


r 

u im 


2  im  ^ jn  ^ in  ^ im ) 


m  jm  t 


(3.25) 


Separation  of  dilatational  and  distortional  components.  In  order  to  disentangle  the 
dependences  of  equations  on  the  Lame  coefficients  A  and  p,  we  separate  the  stress  tensors 
into  its  volumetric  (dilatational)  part,  related  to  the  pressure,  and  distortional  (deviatoric) 
part,  related  to  shear  deformations.  We  decompose  the  stress  tensor  r  as 


Tij  $ij  P  +  aij  i 

(3.26) 

where 

p-.=  -itrr  =  -\rkk 

(3.27) 

is  the  pressure.  The  tensor  a  is  then,  by  construction,  symmetric  and  traceless. 

By  substituting  the  decomposition  (3.26)  in  Eqs.  (3.17)  we  find 

lk  7T  v*  ~  T  dip+  T  di  ay  =  0  ’ 

Po  A0  A0 

(3.28a) 

i  kp  -  (A  +  | p)  dk  vk  =  0  , 

(3.28b) 

i  k  +  p  {di  Vj  +  dj  vt~l  6tj  dk  vk)  =  0  ; 

(3.28c) 

Eq. (3.28b)  is  obtained  by  taking  the  trace  of  Eq.(3.17b).  The  desired  form  of  the  differential 

equations  is  then  simply 

1  k~vi  -  -TdiP+  ~Tdj  aij  =  0  > 

Po  Ao  A0 

(3.29a) 

i  k  pp  -  A0  dk  vk  =  0  , 

(3.29b) 

i  k  +  p  (dt  Vj  +  dj  Vi-l  Stj  dk  vk)  =  0  , 

(3.29c) 
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with  a  dimensionless  parameter 


<P  = 


A  +  -|/i 


(3.30) 


We  did  not  eliminate  the  Lame  coefficient  fi  multiplying  the  derivatives  in  Eq. (3.29c);  had 
we  divided  that  equation  by  /jl,  the  background  medium  limit  /r  — »  0  would  not  exist. 


Equations  in  the  L-S  form.  In  the  following  we  will  represent  Eqs.  (3.29)  in  the  matrix 
operator  form  as 

K,F  =  Q,  (3.31) 

where  the  solution  (field)  vector  X"  is  defined  as 


T  = 


a 


(3.32) 


and  the  action  of  the  operator  1C,  say,  X7  =  ICJ7,  is  represented  as 


4 

*7 

Km 

p' 

= 

cr 

1C  pp 

aij 

XTG\ 

,K/ijk 

Kifki 

(3.33) 


The  above  equation  illustrates  our  indexing  conventions:  elements  of  a  set  of  functions 
( v,p,a )  are  assigned  three  possible  sets  of  indices:  (i,  ,ij)  or  ( k ,  ,kl )  or  (m,  ,mn). 

In  order  to  derive  the  L-S  equations,  we  split  the  operator  1C  as  1C  =  T>  +  V,  where 
V  describes  the  fields  in  the  background  medium,  and  V  the  remaining  interactions.  We 
then  define  the  Green  function  Q  as  the  negative  of  the  inverse  of  the  operator  V,  or,  more 
physically,  as  a  field  generated  by  a  point-like  “unit”  source,  i.e. ,  a  solution  of  the  equation 

VG(r)  =  -152(y)  ,  (3.34) 


where  X  is  an  appropriately  defined  “unit”  tensor.  In  our  case  we  have  to  keep  in  mind 
that  the  symmetric  traceless  tensor  field  a  has  to  be  generated  by  a  source  satisfying  those 
constraints;  hence  the  tensor  X  must  be  symmetric  and  traceless  in  the  indices  associated 
with  the  field  a. 

In  our  case  the  system  of  equations  (3.34)  for  the  Green  function  takes  the  form 


'lk6ik 

-dr 

^ ik 

/->VV 

Vkm 

s? 

r>\G 

y  kmn 

~dk 

ik 

0 

Gpv 

^  m 

qpp 

QP° 

mn 

0 

0 

-se 

3* 

Gap 

ykl 

ncrcr 
y  klmn 

c 

uim 

0 

0 

0 

1 

0 

0 

0 

^ ijmn 

(3.35) 
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where  the  tensor 


^ ijmn  '  2  im  ^ jn  4"  ^ in  ^ jm  3  ^ ij  ^ mn )  (3.36) 

defines  a  symmetric  and  traceless  source  of  the  field  a.  It  is  also  an  idempotent  (or  projec¬ 
tion)  operator, 

^■ijkl  ^ klmn  ~  ^ ijmn  ’  (3.37) 

which  follows  immediately  from  the  tracelessness  property, 


^ iikl  ®  ^ ijkk 


(3.38) 


and  implies 

l'2  =  1.  (3.39) 

Having  specified  the  background  medium  operator  V  (Eq.(3.35)),  we  define  the  inter¬ 
action  term  as  V  :=  1C  —  V.  From  the  system  of  equations  (3.29),  with  Eq. (3.29c)  written 
as 

i  k  atj  +  Aijlk  vk  =  0  ,  (3.40) 

we  find 


yvp 

i 

V3S 

i  k(p-l)6ik 

0 

0 

V  = 

vr 

ypp 

Ki 

= 

0 

i  k  {ip  —  1) 

0 

vap 

*? 

> 

2^  ^ ijik 

0 

0 

(3.41) 


In  order  to  simplify  the  notation,  we  have  temporarily  set  p0  =  A0  =  1;  these  parameters 
will  be  restored  in  the  final  form  of  the  integral  equations. 

We  note  that  the  coupling  between  the  “acoustic  fields”  (v  and  p )  and  the  “shear  field” 
( 7  is  given  here  by  a  differential  operator  proportional  to  /i.  It  is  also  important  to  note 
that  the  Green  function  Q  and  the  operator  (3.41)  satisfy,  by  construction,  the  projection 
conditions 

gi=g=ig  (3.42) 

and 

V J  =  V  =  JV  .  (3.43) 


3.4  An  alternative  form  of  differential  equations  (for  isotropic  media)  in 
the  L-S  form 

We  derive  here  yet  another  form  of  differential  equations,  which  may  be  well  suited  as 
a  basis  of  the  L-S  equation  for  high  contrast  problems.  As  a  guideline  we  will  use  the 
behavior  of  various  physical  quantities  at  material  interfaces  at  which  the  density  (and 
other  parameters)  are  discontinuous:  some  of  the  quantities  are  continuous  across  such 
interfaces,  and  some  not. 
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We  start  with  the  Lame  equation  in  the  first-order,  Eq.(3.8),  and  temporarily  set  p0  = 
= 

k2  p  u,i  +  dj  TtJ  =  0  ,  (3.44a) 

Tij  ~  X  6ij  dk  uk  ~  V  (' di  uj  +  dj  uz)=°  ■  (3-44b) 

On  this  basis,  we  will  first  obtain,  in  a  very  elementary  way,  two  complementary  equations, 
one  for  the  pressure  and  other  for  the  displacement,  in  both  of  which  the  discontinuity  of 
material  parameters  is  “isolated”  in  a  similar  way.  We  will  then  try  to  generalize  those 
equations  to  the  case  of  elasticity. 

Acoustics.  For  acoustics  (p  =  0),  Eqs.  (3.44)  become 

k2  put  +  dj  t1}  =  0  ,  (3.45a) 

Tij~  XSijdiui  =  °  ■  (3.45b) 

By  representing  the  stress  tensor  in  terms  of  the  pressure  p, 

Tij  =  ~SijP ’  (3‘46) 

one  obtains  a  system  of  coupled  first-order  equations 

k2pui-dip  =  0,  (3.47  a) 

p  +  A<9;  =  0  .  (3.47b) 

Acoustics:  L-S  form  of  equations  for  the  pressure.  After  dividing  Eq.(3.47a)  by  p , 
taking  its  divergence,  and  eliminating  the  displacement,  we  obtain  the  usual  second-order 
equation  for  pressure, 

k2p+  \di  Q  dip\  =  0  ,  (3.48) 

where  the  expression  in  the  parentheses  is  continuous.  An  equivalent  L-S  form  of  the  last 
equation  is 

idA  +  fc2)  Q  p)  -  k<2  ( “  “  x)  P  ~  di  (di  ~p)  P  =  °  '  (3'49) 

We  note  two  important  features  of  this  equation: 

1.  The  equation  is  not  really  for  the  pressure  p  (which  is  continuous)  but  for  p/p  (which 
is  discontinuous). 

2.  The  potentially  singular  gradient  of  1/p  is  multiplied  by  p  (which  is  continuous). 
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Acoustics:  L-S  form  of  equations  for  the  displacement.  By  taking  the  gradient  of 
Eq. (3.47b)  and  eliminating  the  pressure  p,  we  find  an  alternative  second-order  equation  for 
the  displacement,  equivalent  to  the  usual  Lame  equation  with  p  =  0, 


k2ui  +  ^diiXdlul)  =  °  • 


(3.50) 


Here,  again,  the  expression  in  the  parentheses  is  continuous.  An  equivalent  L-S  form  of  the 
last  equation,  analogous  to  Eq.(3.49),  is  then 


ididj  +  sij  k2)  (A  uj)  ~  k2  (A  -  p)  ut  -  di  [(dj  A)  Uj\  =  0  .  (3.51) 

Evidently,  Eq.(3.51)  is  (apart  from  the  different  tensor  structure)  analogous  to  Eq.(3.49), 
under  the  replacements 

P  — ►  \  ,  A  — >  -  .  (3.52) 

A  p 

Eq.(3.51)  has  now  properties  analogous  to  those  of  of  Eq.(3.49): 

1.  The  unknown  in  the  equation  is  not  for  the  displacement  ut  (which  is  continuous 
across  interfaces)  but  rather  A ui  (which  is  not  continuous). 


2.  The  possibly  singular  gradient  of  A  is  multiplied  by  Uj,  which  is  continuous. 


Elasticity.  There  is,  of  course,  no  exact  counterpart  of  Eq.(3.51)  in  the  case  of  p  ^  0. 
However,  a  closely  analogous  equation  is 

(' didj  +  5ij  k2)  [(A  +  2/x)  Uj\  -  k2  (A  +  2p  -  p)  Ui  -  di  [(dj  A)  Uj] 

-  dj  [2  di  (p Uj)  —  p  (( di  Uj  +  dj  u-)]  =  0  . 

The  first  three  terms  in  this  expression  have  the  same  form  as  the  corresponding  terms  in 
Eq.(3.51),  and  the  Lame  coefficient  p,  appears  only  in  the  last  term. 


4  Green  functions  in  elasticity 

As  an  initial  step  in  formulating  surface  (boundary)  integral  equations  we  discuss  here 
the  Green  functions  of  the  Larnee  equation  for  the  displacement  field  in  a  general  elastic 
medium.  This  form  of  the  Green  function  will  be  used  in  the  surface  integral  equations 
(Sec.  5). 

On  the  other  hand,  in  the  Lippmann-Schwinger  integral  equations  we  will  only  need  the 
Green  function  for  the  background  medium  -  in  our  case,  air. 

For  completeness,  we  give  here  a  short  derivation  of  the  formulae  for  the  Green  function. 
We  then  represent  it  in  the  form  most  suitable  for  discretization  of  the  integral  equations, 
avoiding  potential  difficulties  associated  with  its  singular  short-distance  behavior. 
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The  Green  function  of  the  Lame  equation.  The  Green  function  of  Eq.(3.1)  -  a 
second-rank  tensor  -  is  defined  by  the  equation 


to2  p  5, 


+  d.n 


( Cimkj  dk )  Gjn(r )  = 


~Sin  ^  (r) 


(4.1) 


and  the  radiation  boundary  conditions  at  infinity.  Physically,  the  function  G-n{ r)  is  the 
displacement  in  the  direction  j  generated  by  a  point-like  force  located  at  the  origin  and 
acting  in  the  direction  n. 

In  the  case  of  an  infinite  homogeneous  medium  the  explicit  form  of  the  Green  function 
(4.1)  is  known  [2,  3],  and  can  be  derived  as  follows: 

We  represent  the  Green  function  G  in  terms  of  its  Fourier  transform  G, 


d  3q 

(27 r)3 


eiqrG„„(q) 


jn\ 


d3q 

(27 r)3 


3iqr 


[/0(?)  6jn  +  / 2(9)  Qj  Qr 


(4.2) 


where  the  Ansatz  for  G  is  based  on  the  assumption  of  the  medium  isotropy.  By  substituting 
Eq.(4.2)  in  Eq.(4.1)  and  solving  for  /0  and  /2  we  find 


Gij  (q)  =  -  ( Su  - 


Qi  Qj\  1  1  Qi  Q 


+ 


1  hj 


t-i  V  13  kl  J  q2  -k |  pkl  q2  -  B 


with 


1 2  2  P 

kr  =  to  - - 

C  A  +  2p 


k2s  =  cu2  P- 
3  P 


The  two  wave-numbers, 

correspond  to  two  sound  speeds, 


k  -  — 

KC  -  _  ’ 

Cc 


ks  - 


UJ 


Cn  — 


'  \  +  2p, 


CS  ~ 


for  longitudinal  (compressional)  and  transverse  (shear)  waves. 
In  the  acoustics  limit  (p  — >  0)  (4.3)  becomes 


1 


GiM)„=0  aa,2  q2  _  ^ 


Qi  Qj 


The  coordinate-space  Green  functions  are  then 

Gb  (r)  =  \  (Sij  +  p  ^  9s(r)  "  ^2  di  dj  9c(r)  , 


and 


^•(r)M=0  =  -  3^2)  iSij  <53(r)  +  di  dj  9C(r))  =  -^pT-p  ( 5ij  6' 3(r)  +  di  dj  9c(r))  > 


(4.3) 

(4.4) 

(4.5) 


(4.6) 


(4.7) 


(4.8a) 


(4.8b) 
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where 


flcW  =  =  ^hoW)  >  5sW  =  4^:  =  ^ho  (V) > 

and  the  spherical  Hankel  functions  of  the  first  kind  (Ref.  [4],  Ch.  10)  are 


h^(z)  =  -~eiz  ,  h^(z)  =  (1  -  iz)eiz  ,  h $\z)  =  (3  -  3iz  -  z2) 

z  zz  z6 


By  using  relations  between  the  Hankel  functions  and  their  derivatives,  we  find 

1  r,  g\r)  ,  {„  g'(r)\\ 

■  +  U  rj  [g  (r)  -  — —  J  j 


—  d ■  d  ■  a(r)  —  —  <f(5- 

^2  i  uj  ;  —  /,2  (/*? 

=  ^  [h51}(fcr)  +  h^(kr)]  +3rirjh ^(Ax)} 


and 


Hence, 


G„(r)  = 


i  kr 


*ij{  1  127T  (A  +  2g) 

i  k c 


h01}(fcCr)  +  (tfy  “  3G  h21)(^Cr) 


127T  g 

In  the  limit  g  — >  0  Eq.(4.13)  becomes 


-2  <*ij  hi1}  (V)  +  ( Sij  -  3  G  f  j)  h21}  (fc 


^(r)  o  = 


i  kc 


^=o  12yrA 
1 


hj)1}  (fccr)  +  (<^-  -  3  r*  fh)  h^  (A;cr) 


Afer  ^ 


Sn  s3(j) , 


(4.9) 


(4.10) 


(4.11) 


5ij  +  ^2  didj)  9(r)  =  ^  {<5y  [2hJ1}(A:r)  -h^(fer)]  +  3  r{  rj  h^Ax)}  .  (4.12) 


(4.13) 


(4.14) 


i.e.,  the  shear- wave  contribution  reduces  to  a  delta-function  term. 

A  dyadic  form  of  the  the  Green  function.  The  Green  function  (4.8a)  or  (4.13)  may 
be  represented  in  a  equivalent  form  involving  dyadic  derivatives.  We  first  express  the  Green 
function  and  its  derivatives  in  terms  of  the  spherical  Hankel  functions  of  the  first  kind, 
(Ref.  [4],  Ch.  10), 


pikr  • 

2(r)  =  (M 

^(r)  =  (M  , 
g"{ r)  =  ( kr )  ; 


(4.15) 
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these  expressions  apply  to  both  the  Green  functions  (4.9).  Hence, 


Vr£>  Vrg(r)  =  1^-  -r®r  g"(r)-^- 

ik2  f  f  hi1*  (At)  „  x  [  (i)"n  .  hi1*  (At)" 
=  —  <  I  -2—4 — -  -  r  (8)  r  hi1'  (At)  -  -2—4 — - 
47t  (hr  [  u  kr 

By  using  relations  between  the  Hankel  functions  and  their  derivatives, 

h<‘>"(2)  =  4‘>W  -  =  hm(z)  +  Ci£i , 

h21}(^)  =  -ho1}(^)  +  ^hi1}(^)  =  -ho1}(^)  -  7ho1}  (z) 

_  47r  r*.2  /  N  ,  oS'(r)  1 


=  -^  ^W  +  3^ 


h'1}(z)  +  hj1}(z)  =  7hS1}(z)  = 


47t  3g'(r) 
i  k2  r 


we  obtain 


h^O)  -  =  ho1}(^)  +  =  ~3Y^  k2g(r)  + 


Vr  (8)  Vrg(r)  =  —  [h°  /fcr)  —  r  <8>  r  h^1} (fcr-)l  , 
47t  At 


/  i  \  47t  T  s  tj'(r) 

r  <g)  r  hi  (A;r)  =  77^  I- - Vr®Vrgr(r)  , 

ifcz  r 


(/  —  3  r  <g>  r)  hi1'*  (kr)  =  hi1*  (kr)  —  3  /  -  —  -  —  Vr  <g>  Vr  g(r) 

lAr  [  r 

=  %  \~k2g(r)  -  6  ^  +  Vr  ®  Vr  5(r) 
i/cz  l  r 

hi1*  (At)  +  (/  —  3r  <8  r)  hi1*  (At)  =  ^  —6  S  ^  +  3  Vr  ®  Vr  ff(r)  , 

ikz  [  r 

—2  hi1*  (kr)  +  (I  —  3  r  <g>  r )  hi1*  (kr)  =  —  3  k2g(r)  —  6  ^  ^  +  3  Vr  <g>  Vr  g(r) 

i  kz  r 


and,  finally, 


G(r)  =  — -2 
A  +  2  g,  y 


(4.16) 


(4.17) 


(4.18) 


+  Vr  (8.  Vr  gc(r)j  -  ^  (-A|  gs(r)  “  2  +  Vr  ®  Vr  </s(r))  . 

(4.19) 


An  alternative  form  of  the  Green  function:  a  reduced  degree  of  singularity.  In 

the  same  dyadic  notation  as  above,  we  can  also  represent  the  Green  function  (4.8a)  in  the 
form 

G(r)  =  I  C(r)  +  Vr  ®  Vr  D(r)  ,  (4.20) 
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with 


C(r)  '■=  77  £'s(7')  >  (4.21a) 

D(r )  ■=  ^2  [#s(r)  ~  9c(r )]  •  (4.21b) 

The  above  representation  exhibits  an  important  property  of  the  Green  function  (4.20): 
since  the  function  D{r)  of  Eq. (4.21b)  is  regular  for  r  — ►  0  (due  to  cancellation  of  the 
singularities  in  the  two  Green  functions),  the  second  term  of  the  Green  function  (4.20) 
is  also  nonsingular,  while,  without  the  cancellation,  it  would  have  contained  a  ~  1/r3 
singularity.  The  reduced  degree  of  singularity  is  particularly  important  in  the  discretization 
of  surface  integral  equations  (Sec.  5). 

5  Surface  integral  equations 

We  discuss  here  briefly  the  form  of  the  surface  boundary  equations  and  their  discretization, 
which  is  now  being  implemented  in  our  solver. 

5.1  Derivation  of  surface  integral  equations  for  scattering  problems 

Surface  integral  equations  in  elasticity  can  be  derived  from  boundary-value  problems  in¬ 
volving  boundary  conditions  defined  on  interfaces  separating  homogeneous  material  regions 
(for  which  unbounded-space  Green  functions,  such  a  discussed  in  Sec.  4,  are  known). 

The  conventional  procedure  is  to  start  with  representation  theorems  expressing  the  fields 
in  a  region  as  an  integral  of  an  appropriate  Green  function  and  the  values  of  the  fields  and 
(possibly)  their  derivatives  on  the  boundary  of  the  region;  typically,  such  representations 
arise  from  the  Green  theorems  and  their  generalizations.  By  imposing  boundary  conditions 
on  the  region-region  interfaces,  one  can  then  obtain  integral  equations  for  the  field  values 
on  those  interfaces. 

More  precisely,  the  above  approach  is  known  as  the  “direct  method”.  Other  forms 
of  integral  equations  can  be  obtained  by  “indirect  methods”  by  postulating  expressions 
(Ansatze)  for  the  fields  in  terms  of  other  sources  supported  on  boundaries  of  the  material 
regions. 

In  elasticity,  a  convenient  form  of  the  representation  theorems  [2,  5]  involves  the  dis¬ 
placement  field  u(r)  defined  in  a  given  domain  fl  and  on  its  boundary  <914,  and  the  traction 
field  f(r)  expressed  in  terms  of  the  stress  tensor  r,  and  defined  on  the  (smooth)  region 
boundary, 

t(r)  :=  n(r)  •  r(r)  ,  (5.1) 

where  n  is  the  exterior  unit  normal  to  dQ. 

In  addition  to  the  second-rank  tensor  Green  function  G  (Eq.(4.20))  for  the  displacement 
field,  the  representation  theorems  also  require  a  third-rank  tensor  Green  function  E  for  the 
stress  tensor;  this  quantity  is  defined  as 

E(r)  :=  A  V  •  G(r)  +  \i  [V  G(r)  +  G( r)  V]  ,  (5.2) 
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or,  in  index  notation, 

V')  :=  AiyS,0,t(: r)  +  ,,  [8,  Gjt(r)  +  dJ  Glt(r)]  ,  (5.3) 

i.e. ,  £  is  related  to  the  displacement  Green  function  G  in  the  same  way  as  the  stress  tensor 
is  related  to  the  displacement  field  (Eq.(3.6)).  We  note  that  the  Green  function  £ is 
symmetric  in  its  first  two  indices. 

It  is  also  convenient  to  introduce  a  second-rank  tensor  Green  function  T  as  a  contraction 
of  £  with  the  normal  vector, 


T(r,  r')  :=  —  h(r')  •  £(r  —  r');  (5-4) 

we  note  that  T(r,r/)  does  not  only  depend  on  the  relative  distance  r  —  r',  but  rather  on  r 
and  r'  separately.  In  analogy  to  r(r,r'),  we  also  write  G(r,  r')  =  G(r  —  r').  Actually,  the 
representation  theorems  involve  the  transpose  Green  functions,  denoted  by  GT  and  Tt. 

It  is  of  interest  to  note  here  that,  since  the  displacement  Green  function  (4.20)  contains 
only  an  ~  1/r  singularity,  the  stress-tensor  Green  function  T  may  contain  at  most  ~  1/r2, 
but  not  ~  1/r3  singularities.  This  fact  facilitates  discretization  of  the  integral  equations 
and  computation  of  the  matrix  elements. 

With  the  above  definitions,  the  basic  form  of  the  representation  theorem,  applicable  to 
a  displacement  field  satisfying  the  homogeneous  Lame  equation  in  the  domain  II  =  I2_,  is 


Ion 


d  V  [r£_  (t,  r')  •  u(r')  +  G]{_  (r,  r)  •  t(r')]  =  < 


u(r)  for  r  G  , 
^u(r)  for  r  G  dQ  , 


(5.5a) 


0 


for  r  G  . 


In  the  second  representation  theorem  the  region  Q_  is  interchanged  with  its  complement 

n+  :=  M3  \  H, 

0  for  r  £  , 

d2^  [rL(t,r')  •  u(r')  +  (r,r')  •  t(r')]  =  -  <J  |u(r)  for  r  G  d£l  ,  (5.5b) 


i  dn 


u(r)  for  r  G  Q  ,  . 


More  precisely,  the  representation  theorem  (5.5b)  assumes  that  the  displacement  field  sat¬ 
isfies  the  Lame  equation  in  and  the  radiation  boundary  conditions  at  infinity;  this  fact 
implies  that  the  displacement  in  Eq.(5.5b)  is  the  scattered  field,  rather  than  the  total  field. 
It  is  also  important  to  remember  that  the  Green  functions  appearing  in  Eqs.  (5.5a)  and 
(5.5b)  correspond  to  different  regions,  =  II  and  £l+.  Finally,  the  expressions  for  r  G  dVL 
have  to  be  interpreted  as  improper  integrals. 

According  to  the  general  “direct  method”  procedure,  the  representation  theorems  (5.5) 
are  now  supplemented  with  the  boundary  conditions  on  an  interface  of  two  solids,  which 
simply  require  continuity  of  the  displacement  and  traction  fields  on  the  boundary  dQ.  It  is 
now  straightforward  to  obtain  the  set  of  integral  equations  for  the  unknown  fields  u  and  t 
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on  the  boundary  9B, 


5  u(r)  - 

1  dV  [r£  +(t,r')  -u(r 

,)  +  G^(r,r').t(r')]=u'"(r)  , 

(5.6a) 

Jdn 

\  u0)  + 

f  dV  [r£_(t,r')  -u(r 
Jdn 

')  +G'iL(r,r')  •t(r')]  =  0  ; 

(5.6b) 

um  is  here  the  incident  field.  Both  equations  hold  at  points  r  £  90  and  in  both  equations 
the  unknowns  are  two  three-dimensional  vectors:  the  displacement  and  the  traction  fields 
u  and  t. 

The  above  form  of  surface  integral  equations  applies  to  the  simple  case  of  a  single  domain 
O  immersed  in  a  background  medium.  It  can  be  generalized  in  a  rather  straightforward  way 
to  the  case  of  multiple  regions  (domains)  separated  by  interfaces.  The  general  form  of  the 
resulting  system  of  equations  is  given  in  Section  2,  Eqs.  (2.1). 

5.2  Basis  functions  and  discretization  of  surface  integral  equations 

In  order  to  solve  the  surface  integral  equations  (5.5)  numerically,  it  is  necessary  to  make 
assumptions  on  the  discretization  of  the  solution,  i.e. ,  on  the  trial  basis  functions,  and  on 
the  test  basis  functions. 

In  our  implementation  we  use  a  discretization  uniquely  determined  by  our  choice  of 
discretization  in  the  volumetric  equations  (discussed  in  detail  in  Sec.  6).  We  also  use, 
similarly  to  the  volumetric  problem,  the  Galerkin  discretization,  i.e.,  identical  trial  and 
testing  basis  functions. 

In  the  volumetric  problem  we  assume  the  displacement  field  is  expanded  in  piecewise 
linear  basis  functions  supported  on  sets  of  tetrahedra.  In  the  corresponding  surface  problem 
we  use,  therefore,  restrictions  of  these  basis  functions  to  the  facets  of  the  tetrahedra  to 
their  boundary  facets  (triangles).  The  resulting  surface  basis  functions  are  piecewise  linear 
vector  basis  functions  describing  the  components  of  the  displacement  field  u.  By  symmetry 
between  the  displacement  and  traction  fields  in  the  integral  equations,  we  assume  analogous 
linear  basis  functions  for  the  components  of  t. 

According  to  the  above  criteria,  we  specify  the  basis  functions  as  follows: 

For  each  vertex  vQ  of  the  surface  mesh  we  define  three  vector  basis  functions,  denoted 
tpa(r),  representing  displacements  in  the  x,  y,  and  z  directions.  Correspondingly,  the  index 
a  refers  to  the  vertex  and  the  direction,  a  =  (va,m),  m  =  1,  2, 3  (or  m  =  x,  y,  z ). 

Each  such  function,  ^(r),  is  associated  with  a  vertex  va  and  supported  on  a  set  of 
triangles  (facets)  fa  sharing  that  vertex.  We  parametrize  the  basis  function  as 

■0a(r)  =  </Va,m(r)  =  emiJr)  ,  (5-7) 

where  em,  is  the  unit  vector  along  the  m-the  axis,  and  <pVa  is  a  scalar  basis  function  defined 

by  _ 

<Mr)  =  Ka(r)  =  ^va,/a(r)  >  (5-8) 

a 

where  the  sum  is  taken  over  the  set  J~a  of  all  facets  fa  sharing  the  vertex  va.  Further,  each 
of  the  linear  functions  4>Va  ja(r),  supported  on  the  facet  fa,  is  uniquely  defined  by  setting 
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its  value  to  unity  at  r  =  va  and  to  zero  at  the  remaining  vertices  of  the  facet.  An  explicit 
expression  is 


^va,/£ 


(r) 


n 


LVa,fa 


Va,  fa 


■  (r 


Xf  (r)  ) 

J  OL  X  ' 


(5.9) 


where  x f,y  (r)  is  the  characteristic  function  of  the  facet  fa ,  hv  j  is  the  unit  outer  normal 
to  the  facet  edge  opposite  the  vertex  va,  and  hVa  ja  is  the  facet  height  relative  to  that  edge. 
Components  of  the  basis  function  ip0.  are  then 


C(r)  =  <a,m(r)  =  <Jmi^Va(r)  •  (5-10) 

As  follows  from  the  construction,  the  scalar  and  vectorial  basis  functions  (5.8)  and  (5.7) 
are  two-dimensional  analogues  (actually,  restrictions)  of  the  piecewise  linear  basis  functions 
supported  on  tetrahedra  and  used  in  the  volumetric  formulation  (Secs.  6.1.1  and  6.1.2).  The 
advantage  of  this  discretization  scheme  is  that  the  solutions  of  the  surface  and  volumetric 
equations  can  be  directly  compared  with  one  another. 


5.3  Structure  of  the  stiffness  matrix  and  computation  of  matrix  elements 

Galerkin  discretization  of  the  integral  equations  (5.6)  results  in  two  types  of  matrix  elements, 


and 


A 


G  _ 

a/3 


[  d2r2^a(ri)  •GT(r1, 


(5.11a) 


(5.11b) 


with  the  Green  functions  corresponding  to  one  of  the  regions  in  question.  The  integrals  are 
taken  here  over  sets  of  facets,  J~a  and  J-  g,  sharing  the  vertices  va  and  v^;  therefore,  they 
can  be  expressed  as  sums  of  integrals  taken  over  pairs  of  facets  fa  and  fp. 

We  give  now  some  examples  of  more  explicit  expressions  for  the  matrix  elements  (5.11b). 
They  are  still  comparatively  simple  relative  to  those  for  the  volumetric  equations,  discussed 
in  detail  in  Sec.  6.  The  procedures  of  simplifying  the  matrix  elements  are  similar  in  both 
cases  and  involve,  mostly,  integration  by  parts  and  using  the  defining  equations  for  the 
Green  functions. 

We  consider  below  the  more  involved  matrix  element  A V ,  in  which  we  set  a  =  ( va ,  m) 
and  /?  =  (vg,  n).  A  contribution  to  this  matrix  element  from  a  pair  of  facets  fa  and  fg  has 
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the  form 


on  f  fi)  —  j/^-va,m;vp,n(.f  on  f (3) 

=  -/  d2^  /  d2r2  <f>  f  (r1)em-rT(r1,r 2)  ■  en4>  (r2) 

jf  jf  va,fa  Vp,Jp 


'fa  Jfp 

]2„  f  i2 


=  -  /  d-Yj  /  d"r2  (j)  (r1)T,inm(r1-r2)ni(r2)(j)  (r2) 

■//a  •//(,  a’U  13,113  (5.12) 

=  -  /  d2rx  f  d2r2(t>Wa  fa{Y1){\8indlGlm{vl-r2) 

J  fa  Jfp 

+  /i  [ch  G 

nm  (r1-r2)  +  a„Gim(r1-r2)]} 

Further  manipulations  of  the  matrix  element  involve  substitution  of  the  representation  of 
the  displacement  Green  function  (4.20),  e.g., 


Glmir  1  -  r2)  =  Slm  °(\rl  ~  r2l)  +  dl  dm.  D(K  ~  r2|)  • 


(5.13) 


One  of  the  resulting  terms,  involving  the  Lame  coefficient  A  and  the  function  D,  becomes 
then 

-\  [  d  2r1  [  d2r2  <f>  (rx)  dt  dt  dm  D(y1  -  r2)  nn(r2)  </>  .  (r2)  .  (5.14) 

_//.  7-  va,/a  V/3-//3 


'fa  Jfp 

This  expression  can  be  further  simplified  by  using  the  relation 

—  Y72 


di  d,  D(r)  =  V2  D(r)  =  gc(r)  -  t  gs(r)  , 


(5.15) 


following  from  Eq. (4.21b)  and  the  defining  equations  for  the  Green  functions  gc  and  gs. 
Finally,  integration  by  parts  transforms  Eq.(5.14)  into 


A  /  /  d2r2»m^„,/„<rl) 


nn(r2)0v  /  fo)  ' 


(5.16) 


Since  the  basis  function  0V  ,  /  is  linear  on  the  facet  fa,  its  gradient  is  the  sum  of  a  constant 
on  that  facet  and  of  linear  delta-functions  supported  on  its  boundaries.  In  the  sum  over 
the  facets  contributing  to  the  full  matrix  element  A^a  the  contributions  of  delta-functions 
from  adjacent  facets  cancel.  This  cancellation  is  complete  in  the  problem  of  a  single  domain 
17,  for  which  dVL  is  a  closed  surface;  however,  delta-function  contributions  from  facet  edges 
may  remain  for  more  complex  topologies  with  several  material  regions. 

Other  terms  in  the  matrix  element  can  be  treated  in  a  similar  way.  For  instance,  the 
term  inolving  A  and  the  function  C{r)  is  very  similar  to  Eq.(5.16), 


^  L  d2ri  If  ^va]  fa  (rl)  Ss  (*T  r2)  ^'n(^"2)  ^Vp,  f  p  (*"2) 


(5.17) 
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The  more  involved  term  proportional  to  the  Lame  coefficient  //  is 


[  d2rJ  d2r  2  di  (/>  (q)  <5mn  ^(q  -  r2)  n,(r2)  ^  (r2) 

■//a  Jfp 

+  [  d\[  d2r 2  dn  <j>  (rl)gc{T1-T2)nm(r2)ct>w  (r2) 

■//«  jf/3  va,/a  VA? 

-75  1  d2rl  [  d2r2diK  f  (rl)dmdn[9c(rl~r2)  -  9s(rl-r2)]ni(r2)<t>v  f(r 2)  ■ 
fcS  Jfa  Jf, 3  “  ^  5 

(5.18) 

Although  the  third  term  in  this  expression  involves  two  derivatives,  the  result  is  well  defined, 
since  the  difference  of  the  Green  functions  appearing  there  is  nonsingular. 


6  Volumetric  integral  equations 

We  describe  here  two  formulations  of  the  volumetric  (L-S)  integral  equations,  based,  respec¬ 
tively,  on  first-  and  second-order  differential  equations  of  elasticity.  In  addition,  we  recast 
the  more  conventional  of  the  equations  into  forms  better  suited  to  handling  high-contrast 
problems. 

On  the  level  of  the  differential  and  integral  equations  themselves,  these  two  approaches 
are  exactly  equivalent.  They  differ,  however,  in  the  choice  of  the  unknowns  and  in  the 
treatment  of  the  material  properties,  hence  in  discretization  aspects. 

In  order  to  allow  direct  comparison  of  the  two  formulations,  we  use  in  both  cases  the 
same  basis  functions;  Therefore,  we  discuss  the  basis  functions  at  the  beginning  of  this 
Section. 

6.1  Basis  functions 

Before  deriving  integral  equations  following  from  the  differential  equations  analyzed  above, 
we  introduce  basis  functions  which  will  be  used  in  the  discretization. 


6.1.1  Piecewise  linear  scalar  basis  functions 


As  an  underlying  form  of  basis  functions  we  will  be  assuming  piecewise  linear  scalar  functions 
supported  on  sets  of  tetrahedra  sharing  a  common  vertex,  and  interpolating  between  1  at 
that  vertex  and  0  at  the  remaining  vertices  of  the  set  of  tetrahedra.  We  represent  such  a 
function  <pa  as 

</>a(r)  =  (r)  =  t  (r)  ’  O3'1) 

V  Ot  <  ^  V  OLt^OL 

taeTa 

where  the  sum  is  taken  over  the  set  T a  of  all  tetrahedra  ta  sharing  the  vertex  va,  and  each 
of  the  linear  functions  4>Va  t  (r),  supported  on  the  tetrahedron  ta,  is  defined  to  be  unity  at 
r  =  va  and  zero  at  the  remaining  vertices  of  ta.  An  explicit  expression  is 


1  - 


n 


y,tQ 


va,t0 


(r-  vc 


XtA  r)  , 


(6.2) 
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where  Xta(r )  is  the  characteristic  function  of  the  tetrahedron  ta,  nv  t  is  the  unit  outer 
normal  to  the  tetrahedron  face  opposite  the  vertex  va ,  and  hVa  t  is  the  tetrahedron  height 
relative  to  that  face.  When  using  such  basis  functions  we  assume,  as  before,  that  the 
material  parameters  are  constant  on  the  tetrahedra. 

First  derivatives  of  the  basis  function  cj)a.  The  basis  function  <pa  is,  obviously,  contin¬ 
uous  and  differentiable  (except  for  interfaces  between  the  tetrahedra,  where  the  derivatives 
do  not  exist).  However,  the  individual  functions  cpVa  ta  are  not  continuous,  because  they 
characteristic  functions  are  not;  therefore,  their  gradients  are  constant  inside  the  tetrahe¬ 
dra,  but  involve  also  delta-functions  concentrated  on  the  interfaces  between  the  tetrahedra 
in  the  set  T a.  Explicitly, 

V<^va,4r)  =  -J^—Ka,ta  Xta  (r)  —  ^vcta  (r)  n(r)  6gta  (r) 

V ot  i  ta 

=  -]^—nva,taXta(r)-  E  ^v«,/a(r)  4  W 

Va’ta  f*edta 

Here,  in  Eq.(6.3a),  59ta  is  the  delta  function  concentrated  on  the  tetrahedron  boundary, 
and  n(r)  is  the  outer  unit  normal  to  that  boundary.  In  Eq.(6.3b),  fa  E  dta  are  faces  of 
the  tetrahedron  ta ,  n^^  is  the  unit  normal  to  the  face  fa  in  the  direction  out  of  the 
tetrahedron  ta ,  and  (i>Vaja  is  the  basis  function  4>Vatta  restricted  to  the  face  fa.  Because  of 
continuity  of  cj)a,  the  function  <^>Vq  ^  is  well  defined;  it  is  linear,  and  interpolates  between  1 
at  the  vertex  va  and  0  at  the  remaining  vertices  of  the  face. 

The  delta-function  contribution  to  the  gradient  (6.3)  vanishes  on  the  exterior  boundary 
of  the  set  T a  of  the  tetrahedra  (because  the  basis  function  <j)Va  ta  vanishes  there),  but  it  is 
nonzero  on  the  other  faces. 

The  delta- function  terms  cancel  pairwise  in  the  sum  over  the  tetrahedra  ta  E  T Q ,  hence 
the  gradient  of  function  (f>Va  is  regular,  consistently  with  the  fact  that  the  function  <f>Va  itself 
(given  by  the  sum  of  (6.1))  is  continuous. 

However,  the  delta  functions  may  give  nonzero  contributions  if  the  vertex  va  is  located 
on  the  object  boundary  dQ,  and  thus  some  of  the  facets  fa  are  not  interfaces  but  rather 
the  boundary  facets  of  the  object. 


(6.3a) 

(6.3b) 


Second  derivatives  of  the  basis  function  <f)a.  Several  expressions  for  matrix  elements 
(e.g.,  (6.59d)  and  (6.59e)  in  the  following)  involve  second  derivatives  of  the  basis  function 
0V  ,  which  are  delta  functions  and  derivatives  of  delta  functions,  concentrated  on  the  in¬ 
terfaces  of  the  tetrahedra  in  the  set  Ta  (i.e. ,  tetrahedra  sharing  the  vertex  va).  From  the 
expression  (6.3b)  for  the  first  derivatives  we  obtain,  for  each  of  the  basis  functions  supported 
on  a  tetrahedron  ta, 


dm  dl  <Ava,ta(r)  =  T -  E  4  +  <>ta  */>) 

Va’ta  fa£dta 

-  E  "£  "/a  ,  fa  (r)  4  (4  ’  r )  • 

fu&dta 


(6.4) 
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The  notation  used  here  for  the  delta  function  derivative  actually  means 


6U  (“/« > r)  =  Xfa  (rfa )  S'fa  (hfa  •  r)  ,  (6.5) 

where  r ^  is  the  projection  of  r  on  the  facet  fa  along  the  direction  .  The  derivatives  of 
delta  function  contribute  only  for  facets  fa  having  va  as  one  of  their  vertices,  since  on  the 
facet  opposite  the  vertex  va  the  function  <t>Vaja( r)  vanishes. 

In  the  double  derivatives  of  the  basis  function  ef>  (r)  (Eq.(6.1))  each  interior  facet 
fa  ^  a  contributes  to  two  adjacent  tetrahedra,  and  the  derivatives  of  the  delta  functions 
cancel,  since  r)  +  d^(— nja,r)  =  0.  For  boundary  facets  (fa  G  dQ),  however,  the 

derivatives  of  the  delta  functions  remain;  hence 


dm9^Va(r)=  J2  T -  E 


«,,*«  +  Ka,ta  n?a)  SfJr) 


taeTa  v“-‘“  fa£dta 
fa  &Fa  n  dVL 

-  E 

fa^%a 


1 


h 


Va ,  ta  + 

1 


h 


«,,*»+"/«  +"ta,ta+«/L) 


Vq,,  tc- 

E  "/«  fa  (r)  4  (“/a  >  r)  • 


/a6/'an9S1 


(6.6a) 


(6.6b) 


The  second  form  of  this  expression  is  obtained  by  grouping  contributions  from  each  facet 
fa  and  denoting  by  ta+  and  ta—  tetrahedra  on  the  positive  and  negative  sides  of  the  facet 
(according  to  the  direction  of  the  normal  ) .  The  sum  in  Eq.(6.6b),  however,  is  taken  over 
all  facets  fa  belonging  to  any  tetrahedron  in  the  set  Ta,  not  only  over  facets  fa  G  T a 
sharing  the  vertex  vQ;  we  indicate  this  by  writing  fa  G  Ta.  If  the  facet  fa  is  an  exterior 
facet  of  the  set  Ta,  i.e. ,  fa  G  dTa,  the  contribution  of  the  tetrahedron  ta+  is  absent. 

A  simple  geometrical  analysis  shows  that,  for  interfaces  fa,  the  linear  combination  of 
the  normals  nVa,ia+  and  in  Eq.(6.6b)  is  proportional  to  hjQ , 


nv  +  _i_  nv  +  _ 

v  ot  i  i'Ot  i  v  a  i  lol 


~  n. 


h  h  *  /« 

LVa,ta+  lVa,ta- 


(6.7) 


Therefore,  Eq.(6.6b)  can  be  written  as 


dmd^vJr)=-2  E 


n  f  ,  • nf  n  f  _ ■ nf 

v  ot  i  i  Joe  v  at  >  J  a 


K 


f  (zTf  \  va,tQ  + 

J  a  v=./  ex. 


Lva,ta- 


hT  h\  5  f  (r) 

Joe  Jol  Jax  ' 


+  2  E 


1 


_  h  hTa^fa6fa^ 

fad&Ta  Va’ta 

E  ™Ta  "/a  6fa  (“/„  >  r)  Ka,  fa  (F)  ' 

faerandsi 


(6.8) 


We  have  separated  here  contributions  from  the  facets  sharing  the  vertex  va  and  the  facets 
on  the  boundary  dTa  of  the  set  of  tetrahedra.  In  the  latter  contribution  hVa  t  is  the  height 
of  the  tetrahedron  adjacent  to  the  facet  fa,  measured  from  the  vertex  va. 
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6.1.2  Piecewise  linear  vector  basis  functions 


In  the  following  we  consider  linear  vector- valued  basis  functions  i/?a(r)  associated  with 
vertices  and  supported  on  sets  of  tetrahedra  sharing  the  given  vertex.  For  each  vertex 
va  we  define  three  vector  basis  functions  representing  displacements  in  the  x,  y,  and  2 
directions;  the  index  a  in  if)a  refers  thus  to  the  vertex  and  the  direction,  a  =  (va,m), 
rn  =  1,2,3  (i.e.,  m  =  x,  y,  z ). 

We  parametrize  these  functions  as 

^a(T)  =  </Va,m(r)  =  em  ijr)  ,  (6-9) 

where  em  is  the  unit  vector  along  the  m-the  axis,  and  (f>Va  is  the  scalar  basis  function  defined 
by  Eq.(6.1).  Components  of  the  basis  function  (as  appearing  in  Eqs.  (6.58))  are  then 

C(r)  =  <0,m(r)  =  5 mi  0v„(r)  •  (6-10) 

The  derivatives  of  the  basis  functions  appearing  in  Eqs.  (6.58)  represent  either  the 
pressure  or  the  strain  tensor, 

*va,m(r)  =  3>va,m(r)  =  dmKSr )  (6-n) 

or 

®v0lm(r)  =  5  [9l  ^vQ,m(r)  +  9j  <a,m(r)] 

=  5  6 rni  93  <t>va  (r)  +  5mj  dl  (j)Va  (r)]  ;  ^  ^ 

these  expressions  are  valid,  of  course,  for  any  (not  only  linear)  scalar  basis  functions  < f>.  We 
note  that  the  divergence  of  the  basis  function  representing  displacement  in  the  direction  m 
is,  actually,  a  derivative  of  the  scalar  basis  function  in  that  direction. 

6.2  Integral  equations  in  first-order  formulation 

General  structure  of  the  L-S  equations.  The  system  of  equations  (3.29)  with  a  pre¬ 
scribed  source  S  has  the  form 


KT=  (V  +  V)F=  -S  ,  (6.13) 

where  the  source  is  assumed  to  act  in  the  space  of  traceless  symmetric  tensor  fields  a, 

IS  =  S  .  (6.14) 

It  can  be  now  verified,  with  the  use  of  projection  properties  (3.39),  (6.14),  (3.42),  and  (3.43), 
as  well  as  the  definition  (3.34)  of  the  Green  function,  that  Eq.(6.13)  can  be  represented  as 
the  L-S  equation 

{l-gV)F  =  GS=:Fh\  (6.15) 

where  is  the  incident  field. 
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It  follows  immediately  from  Eq.(3.35)  that  some  blocks  of  the  Green  function  are  iden¬ 
tically  zero, 


gvv 

Gvp 

gva 

Gpv 

qpp 

gp* 

0 

0 

gaa 

(6.16) 


Therefore,  the  kernel  Q  V  appearing  in  the  L-S  equation  (6.15), 


gvv 

gvp 

gva 

Vvv 

0 

0 

yyvv 

yyvp 

0 

W  :=GV  = 

gpv 

gPP 

gp« 

0 

vpp 

0 

= 

yypv 

Wpp 

0 

0 

0 

gva 

yav 

0 

0 

wav 

0 

0 

also  exhibits  a  number  of  zero  blocks. 


(6.17) 


Evaluation  of  the  Green  function.  The  system  of  equations  (3.35)  for  the  Green 
function  components  can  be  easily  solved  in  terms  of  the  Fourier  transforms 

GabXr)  =  J  7^j3  eiq'r<?“fc(q)  ,  o,6  =  v,p,a.  (6.18) 


Eq.(3.35)  becomes  then 


[k6ik 

-i  Qi 

i5ik<h 

“i  % 

i  k 

0 

0 

0 

ikSikSjl 

~VV 

y  km 

g7 

~va 

y  kmn 

?>pv 
y  m 

^p 

GV<T 

mn 

(q) 

~crv 

~crp 

~<J(T 

g klm 

g m 

y  klmn 

(6.19) 


c 

uim 

0 

0 

0 

1 

0 

0 

0 

ijmn 

By  considering  the  last  row  of  equations  in  the  system  we  find  immediately 

~<TV 

^  klm^Si)  0  , 

Slf(q)  =  0  , 

~<TCT  .  .  1 

GklmnvX)  ~  ^klmn  ' 


(6.20a) 

(6.20b) 

(6.20c) 
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After  substituting  these  expressions  in  the  remaining  equations,  we  obtain 


i  k  0Zi( q)  -  i  Qi  &m  (q)  =  -  dim  >  (6.21a) 

i  k  aJP(q)  -  i  £PP(q)  =  0  ,  (6.21b) 

~V(T  .  .  ^XPCT  .  .  1  . 

1  k  Gimn{ q)  -  1  qt  Gmn{ q)  +  1  Ql  Jj  Ailmn  =  0  >  (6-21c) 

- — '  VV  - — 'DV 

-  i  Qk  Gkm( q)  +  i  k  Gm  (q)  =  0  ,  (6.21d) 

-i<?fc<?fcP(q)  +  ifc^PP(q)  =  -  1  ,  (6.21e) 

-  i  Qk  GZnni q)  +  i  k  £q„n(q)  =  0  •  (6.21f) 


The  subsystem  of  equations  (6.21a),  (6.21b),  (6.21d),  and  (6.21e)  is  associated  with  the 
acoustic  problem,  described  by  the  fields  v  and  p.  It  can  be  easily  solved  by  substituting 
the  Ansatze  based  on  rotational  invariance, 

—  VV 

Gkmi q)  =  a(g)  8km  +  6(<?)  Qk  Qrn  ,  (6.22a) 

^IP(q)  =  GVk  { q)  =  c(q)  qk  ,  (6.22b) 

GPP{ q)  =  d(q)  .  (6.22c) 


The  solution  is  then 


G  kmG l) 
GlP(  q) 
^P(q) 


i 

k 


_Qk_Qm_\ 
q2  —  k2  ) 


epv(q) 


~i  Qk 
q2  —  k2 


— i  k 
q2  —  k 2 


(6.23a) 

(6.23b) 

(6.23c) 


Out  of  the  remaining  equations,  Eq.(6.21f)  allows  us  to  express  QP  in  terms  of  G 


/  \  Qi  p,v<T  /  x 
Gmn\ q)  —  Glmn{ q)  , 


and  the  last  equation  to  be  solved,  (6.21c),  becomes 
(■ k 2  &u  -QiQl)  0lmn( q)  =  -  i  Ql  Ailmn 


—  2  (^im  Qn  +  8 in  Qm  3  Qi  &mn) 


(6.24) 


(6.25) 


By  substituting  here  the  Ansatz 

GZn( q)  =  ^(?)  Aklmn  %  +  Aijmn  ?i  Qj  Ql 

— vcr  — per 

we  find,  finally,  solutions  for  Q  and  then  G  as 

'P,Yry  /x  1  (  *  Aijmn  QiQ  j  Ql  \ 

Glmn\ q)  —  ~~  ^2  yAklmnQk  q2  _  £.2  J 

piP*7  /x  1  Aijmn  Qi  Q  j 
5'””<q)  =  k  q2-fc2  • 


(6.26) 


(6.27a) 

(6.27b) 


28 


These  expressions  are  manifestly  symmetric  and  traceless  in  the  indices  m,  n. 

The  structure  of  the  Green  function  matrix  and  its  blocks  (Eqs.  (6.23),  (6.20c),  and 
(6.27))  can  be  summarized  as 


£(q) 
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q7 

~va 

y  kmn 

5pv 

qpp 

GPa 

y  mn 
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y  klm 
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yki 
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y  klmn 

(q) 


i  (sr  QkQm  \ 

-i  Qk 

i  (  a  ^ ijmn  Qi  Qj  Qk  \ 

k  \  km  q2  —  k2  ) 

q2  —  k 2 

f, 2  l  klmn  ^1  q2  _  £.2  J 

-i  Qm 

— i  k 

i  ^ ijmn  Qi  Qj 

q2  —  k 2 

q2  —  k 2 

k  q2  —  k2 

0 

0 

kAklmn 

(6.28) 


The  Green  function  in  coordinate  space.  Fourier  transformation  (6.18)  of  the  Green 
function  operator  (3.41)  back  to  the  coordinate  space  involves  just  the  substitutions  qk  — > 
— i  dk  and  the  usual  Helmholtz-equation  Green  function 


g{ r)  =  g(r)  ■■= 


d3c/ 


j  q  r . 


1 


Jikr 


(27r)3  q2  —  k2  —  iO  4nr 

(with  r  :=  |r|),  satisfying  the  equation 

(V2  +  k2)g(  r)  =  —  <53(r) 


(6.29) 

(6.30) 


and  the  radiation  boundary  condition.  Hence, 
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y  kmn 
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qpp 
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mn 

nc rv 
^klm 
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yki 
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y  klmn 

l  [4m  63 (r) 

+  dkdm.g(  r)] 

-  dk  44 

“P  [AklmndlS3(r) 

+  Aijmndidjdkg(r )] 

-  dm  g{ r) 

-i  kg(r) 

-  1  Aijmn  3idj  44 
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lAklmJ3^ ) 

By  using  the  definitions  (3.36)  and  (6.30)  of  the  tensor  A  and  the  Green  function  g ,  one 
can  further  simplify  the  blocks  Qva  and  Qpa  to 


5L(f)  =  -  k  2  [\  (4m  3n  +  4n  3m)  4(f)  +  4  ffmnW]  -  (6.32a) 

5»(r)  =  -  i  k~l  (3  5m.n  4(r)  +  gmn )  ,  (6.32b) 
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where  we  defined,  for  convenience,  an  auxiliary  symmetric  Green  function 


9mn(r)  ■=  (|  Smn  +  dmdn)  S(r)  ■  (6-33) 

The  L-S  integral  equations.  An  explicit  form  of  the  integral  equations  (reflecting  the 
structure  of  Eq.(6.17))  is 

P( r)  Vi(r)  +  J  dV  [didmg{ r  -  r')]  [p(r')  -  l]  vm(r') 

+  ^2  j  d3r-/  {  [ d m  53(r  -  r')]  M(r')  [d'i  vm(r')  +  d’m  v^r’)] 

+  [di9mn(r  ~  r')]  2/x(r')  d'mvn(r')} 

+  i k  J  d3/  \di  g( r  —  r')]  [<p( r')  —  l]  p( r') 

=  ^n(r)  ,  (6.34a) 

p( r)  —  k2  J  d3r'  g(r  —  r')  \p(r')  —  l]  p( r') 

+  -^  p{r)dmVm{r)  +  j_  J  d3r' gmn(r  -  r')2p(r')  d'mvn(r') 

=  pin(r)  ,  (6.34b) 

<Tij(T)  -  j:  Mr)  [di  Vj(r )  +  dj  Vi(r )  -  |  Sij  °m  Vm(T)] 

=  Vij(r)  ■  (6.34c) 

In  particular,  because  of  the  vanishing  third  column  of  the  kernel  VV  of  Eq.(6.17),  the  stress 
tensor  a  does  not  appear  at  all  in  the  first  two  equations,  (6.34a)  and  (6.34b).  These 
equations  can  be  thus  solved  for  v  and  p.  and,  if  desired,  the  stress  tensor  can  be  evaluated 
by  using  Eq.(6.34c),  which  just  reproduces  Eq.(3.29c). 

The  integral  equations  (6.34)  above  are  written  in  the  form  assuming  p0  =  A0  =  1.  After 
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reintroducing  these  parameters,  equations  (6.34)  take  the  final  form 


P(r) 

Po 


;i(r)  +  j  dV  (didmg{ r  -  r'))  -  1^)  vm(r') 


+  ¥dr‘ 


M(r) 


idivmir)  +  dmvi(r)) 


Mr')  a/ 


+  p  J  d3r'  (di  9mnir  ~  r0)  9m  Vn(r') 

+  Y  J  dV  (a.  ff(r  -  r'))  (<M)  -  l)  p(r') 

=  , 

p( r)  —  k2  f  d3r'  g(r  —  r')  (</?(r')  —  l)  p(r') 


+ 


2i  At(r) 
3/c  A0 

=  Pin(r)  , 


um(r)  +  -  dr  gmn( r  -  r  ) 


^  Mr')  a, 


An 


f  r / ) 


%  (r)  -  l  ^  ft  ^(r)  +  dj  vi(r)  ~  |  ^  °m  Mr)] 
=  o§(r)  . 


(6.35a) 

(6.35b) 

(6.35c) 


6.2.1  Matrix  elements:  general  expressions 


We  discretize  the  L-S  equations  (6.35)  by  using  the  Galerkin  method,  in  terms  of  vector 
and  scalar  basis  functions  rj)a  and  <f>a  for  the  fields  v  and  p.  The  resulting  stiffness  matrix 
takes  then  the  form 


A  = 


Avv  i  ylvp 
ylpv  |  App 


(6.36) 
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with  the  blocks 


Al}=  [  #r*(r)h£)^(r) 

J  P  0 

-  J  d 3r1  j  d3r2  (<9!C(r,))  (<9f'  g(r,  -  r2))  “  l)  V$(r2) 

-  p  /  d3r  (0^(r))  ^  (0^(r)  +  d^(r)) 

-  p  y  d3ri  J  d3r2  (31C(ri))  £fc/(ri  -  r2)  (^2^(r2))  ,  (6.37a) 

A7p  =  -yo  j  d3^i  J  d3r2  (d!Y>a(ri))  g(rx  -  r2)  (p(r2)  -  l)  <^(r2)  ,  (6.37b) 

A^  =  ^k  /d3r<^«(r)M(r)  (^'^(r))  (6.37c) 

+  ^T  /d3ri  f  d3r2(t)a(ri)9ki(r1  -r2)n(r2)  (dfylp(r2))  ,  (6.37d) 

=  ~k2  J  d3ri  J  d3r2  4*(ri)  ff(ri  -  r2)  Mr2)  -  1)  </>/3(r2)  •  (6.37e) 

6.2.2  Matrix  elements  with  composite  linear  basis  functions 


In  the  following  we  assume  piecewise  linear  vector  and  scalar  basis  functions,  as  described 
in  Sections  6.1.2  and  6.1.1.  By  expressing  the  vector  basis  functions  in  terms  of  the  scalar 
ones  (Eq.(6.9)),  we  can  represent  Eqs.  (6.37)  as 

A7a,m-^p,n  =  Smn  \  d3r  —  Ka(r) 

J  P  0 

-  J  d3^1  J  d3r2  (d?  g(r1  -  r2))  -  -  1^  </>V/3(r2) 

“^2  /d3r^  Smn{&Ka{r))  +  {dnKa(r)) 

~  ^2  j  d3ri  f  d3r2  (^v.^i))  9kn(r  i  -  r2)  (92</>v5(r2))  ,  (6.38a) 

Ava,,m;vp  =  ~yo  f  d3ri  /  d3r2  (9r^va(ri))  #(ri  -  r2)  (^(r2)  -  l)  </>V/3(r2)  ,  (6.38b) 

=  §:  J  d3r  cj)Va(r)  g(r)  (3>V/3(r)) 

+  jf  J  d3^2^vQ(rl)fffcn(rl  -r2)^(r2)  (520v5(r2))  >  (6.38c) 

^v!(;v3  =  ~k2  f  d3ri  J  d3r2  ^va  (ri)  fl,(ri  -r2)  (</?(r2)  -  1)  0V/S(r2)  5  (6.38d) 

as  before,  vector  basis  functions  are  labeled  by  pairs  of  indices  (vQ,  m)  or  (v^,  n),  represent¬ 
ing  the  vertex  and  the  Cartesian  component  of  the  vector,  and  the  scalar  basis  functions 
are  indexed  by  the  vertices  only. 
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6.2.3  Matrix  elements  with  elementary  linear  basis  functions 

As  before,  we  now  have  to  express,  in  Eqs.  (6.38),  the  “composite”  basis  functions  (sup¬ 
ported  on  sets  of  tetrahedra)  in  terms  of  “elementary”  basis  functions  supported  on  individ¬ 
ual  tetrahedra  or  facets.  By  using  the  relations  of  Section  6.1.1  we  find  the  displacement- 
displacement  matrix  elements  in  the  form 
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(6.39a) 


(6.39b) 


(6.39c) 


(6.39d) 


(6.39e) 


(6.39f) 


where  va  is  the  volume  of  the  tetrahedron  ta.  The  last  two  terms,  Eqs.  (6.39e)  and  (6.39f), 
can  be  integrated  by  parts  and  rewritten  in  a  form  involving  only  first  derivatives  of  the 
Green  function,  e.g., 


/ 


(6.40) 


After  substituting  this  identity  in  Eqs.  (6.39)  we  would  find  that  some  contributions  of  the 
facets  fp  cancel:  this  happens  whenever  the  two  tetrahedra  t  adjacent  to  the  face  have  the 
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same  Lame  coefficient  //(f).  In  other  words,  the  contributions  to  the  r2-integral  come  only 
from  the  discontinuities  of  //.  We  can  obtain  the  same  result  by  integration-by-parts  in 
Eq. (6.38a),  which  would  pick  derivatives  of  //,  hence  delta-function  contributions. 

Similarly,  the  displacement-pressure  and  pressure-displacement  matrix  elements  are 


^;V/3  =  V-  E  E  MV  -  !) 


toL^'Ta 


K  t  Wq 

v  a  5  t'cx. 


f  d3rj  d 3r2^^(r2)S(rrr2) 

J  t(x  w  t/Q 


^  fa&FandQtpGTf} 


/  d\  d3r2  Va,/Q(ri)  -r2) 

'/a  Jtp 


(6.41a) 


(6.41b) 


4pv  =  — 

va  ;vg,n 


(§  EZ  /  d3r<(>Wa(r) 

ic«eTa  v/3da 

I  E  E 


tadTa  tptz'Tp 


(6.42a) 


[  d3r1  /  d3r2</>Va >ta(r1)  (dfd?  ^  -  r2))  .  (6.42b) 

W  t(X  **  t p 

The  second  derivatives  of  the  Green  function  can  be  eliminated,  as  in  Eqs.  (6.39e)  and 
(6.39f),  by  using  the  identity  (6.40). 

Finally,  the  pressure-pressure  matrix  elements  are  given  by 


AZn  =  -k2  £  £  (*>(*„) - 1) 

toL^-Ta,  tpeTp 


/  d3ri  /  d3r2  <t>Vatta  (ri)  </>V/3,t/3(r2)  <?(iT  —  r2) 

«y  ^  Q.  i/ 


(6.43) 


6.2.4  Summary  of  the  expressions  for  the  “basic”  matrix  elements 

The  basic  matrix  elements  appearing  in  Eqs.  (6.39)  -  (6.43)  are  as  follows: 

tetrahedron-tetrahedron  matrix  elements: 
constant-linear: 

TT1  A(ta;vp,  tg)  =  I  d3r1  I  d3r2  (j)Vp  tp{r2)  g(r1  -  r2) 

J  tcx  **  t  p 


TT2  A\ta\\p,  tg)  =  /  dq  d  r2  <f>  t  (r2)  dl  g(r1  -  r2) 
Jta  Jtp 
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linear-linear: 


TT3  A(va,  ta;vg,  t0)  =  J  d3r1  J  d3r2  ^Va,ta(r1)  0V/8>t/S(r2)  gfa  -  r2) 


Up 

tetrahedron-facet  matrix  elements: 
constant-constant: 

TF1  A\ta;fg)=  f  d3r1  [  d2r2  dl  g(r1  -  r2) 


linear-constant: 


TF2  A'(va,ta-,f0)=  j  d3r±  j  d2r2  ^  ^(rj  d^(r1  -  r2) 


linear-linear: 


TF3  A(va,  tQ;vg,  f0)  =  d3ri  d2r2  ^^fo)  </>V/3i//3(r2)  g{r1  -  r2) 
TF4  A\va,ta]Vg,  fg)=  j  d3rx  f  d2r2  <^Va,ta(r1)  0V/3,//3(r2)  ^  s(iT  -  r2) 


facet-facet  matrix  elements: 
linear-constant: 


FF1  Al(vaJa-Jg)=  d2rx  /  d2r2</.Vaita(r1)c)*5(r1  -  r2) 


fp 


Whenever  derivatives  of  the  Green  function  appear,  they  are  not  normal  derivatives. 


6.3  Integral  equations  in  second-order  formulation 

Eqs.  (3.9)  have  the  general  form  of  an  “L-S-type”  differential  equation  with  three  alternative 
expressions  for  the  differential  operator  V.  It  can  be  easily  verified  that  the  background- 
medium  Green  function  g  =  — ZW1  is  given  by 


9ij {r)  =  ~k  Sti  5  (r)  +  di  d-  g{ r) 


*  j  • 


and  satisfies  the  equation 


(■ k 2  Sij  +  di  dj)  gjl(r)  =  -  Su  53( r) 

and  the  outgoing- wave  boundary  conditions.  The  resulting  L-S  equation 


(6.44) 


(6.45) 


«i(r) 


r')Vij(r,)uj(r') 


(6.46) 
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becomes  then,  with  Eq.(6.44) 


Uj(r)  +  fc  2Vij(r)u.j(r')  +  k  2  /  dV  g(r  -  r')  VH(r')  =  'u-n(r/)  .  (6.47) 


Depending  on  the  form  of  the  differential  equation,  various  forms  of  the  L-S  equations  can 
be  obtained.  In  the  following  we  consider  three  such  forms,  corresponding  to  Eqs.  (3.9a), 
(3.9b),  and  (3.9c),  and  we  refer  to  them  as  equations  of  type  (a),  (b),  and  (c). 

The  interaction  operator  V  derived  from  the  differential  equation  (3.9a)  gives  rise  to  an 
L-S  equation  of  the  “basic”  form 


ui(r)  -  (  1  -  —  )  ui( r)~k2d. 
P  o 


1  -  )  d.uAr) 


An 


+  k 2  9, 


M  r) 


diuj(  r)  +  djUi{r)^ 


—  J  d3r'  di  dt  g(r  —  r')  ^1  —  ^  ui{r') 

—  AT2  J  d 3r'  di  dl  g{r  —  r')  d[ 

+  kT2  J  d 3r'  di  dl  g( r  —  r')  dj 
=  MM  • 


1  - 


MO 


An 


djUji O 


d'iuj{  r')  +  d'jUl(  r')) 


(6.48) 


With  the  operator  V  corresponding  to  Eq.(3.9b)  we  obtain 


uAr)  +  k  2  ~TT  dj 

P{  r)  3 


/r(r) 

An 


diUj{r)  +  djUi(r)^ 


-  k  2  f  ^a.wir) 


p(r)J  A0  3  3 


1  - 


PoA(r') 


djuj{ O 

P(r' 


+  j  dVfl^r-rO  f  Aop(r/ 

+  AT2  J  d3r/  9,.  dl  g( r  -  r')  -^y  S' 

-  k~2  j  dV  9.  c)z  c/(r  -  r')  ^  dj  Uj(r') 

=  MM  , 


^^■(O  +  ^«i(r')) 


(6.49) 


where  <9  and  &  denote  derivatives  with  respect  to  r  and  r'. 

In  deriving  Eq.(6.49)  we  applied  integration  by  parts  to  the  term 


-k  2  J  d3/  di  d[  g( r  —  r')  d[ 


resulting  from  directly  from  the  original  L-S  equation  of  (6.47) .  We  apply  here  the  derivative 
di  to  the  Green  function,  and  use  its  defining  equation  to  obtain  the  fourth  term  in  Eq.(6.49) 
with  a  single  derivative  of  the  Green  function. 

The  boundary  term  in  integration  by  parts  vanishes  because  the  material-dependent 
term  (in  the  square  brackets)  vanishes  outside  D. 
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Finally,  with  the  operator  V  of  Eq.(3.9c)  the  L-S  equation  takes  the  form 


+  AT 2  d 


^(r)  [diuj(r)  +  diui(r 


j 


-  k-z  (a  t  ',(r) 


M1-)/  Po 

£A(r)  Sij  dkuk(r )  +  ^(r)  (^(r)  +  dj«i(r 

+  J  d3r'dig{r  -  r')  [l  -  £A(r')]  dkuk(r') 

+  k~2  J  dV  <9- <9,#(r  -  r')  <9'  ^(r')  (^ui(r')  +  dju,(r')) 
-  AT2  J  d V  a.  8t  g{v  -  r ') 

PA(r')  ^9fcufe(r/)  +  ^(r/)  (^(rO  +  d'u^r' 

=  <(r)  , 

or,  after  some  rearrangements  and  integration  by  parts, 

«i(r) 

-k~2dj 


(!  -  ^A(r))  Sijdkuk(T)  -  ^(r)  5*«i(r)  +  ^*( 


-  (  3,  Po 


Pi  r) 


0(1) 

0(P0/P) 

0(1) 


|pA(r)  <5ij  5fc%(r)  +  PM(r)  +  djui{r 

—  k~2  J  d3r'  di  dt  d-  g(r  —  r') 

(!  -  £A(r0)  Stj  d'kuk{ r')  -  ^(r')  (d{Uj(r')  +  d-u,(: r'))  0{pj p) 

-  k~2  j  dV  8i  fJt  <?(r  -  r')  (at 

VXW)  5ij  d'kuk(r')  +  ^(r/)  (d'iuj(r')  +  djui(r' 


=  «F(r)  • 


0(1) 

0(1) 


(6.50a) 

(6.50b) 

(6.50c) 

(6.50d) 

(6.50e) 


(6.50f) 

(6.50g) 

(6.51a) 

(6.51b) 


(6.51c) 


(6.51d) 


(6.51e) 

(6.51f) 


In  the  last  forms  of  the  L-S  equation  we  indicated,  next  to  each  term,  its  expected  magnitude 
in  the  high  contrast  limit  p0/ p  —>  0  (the  origin  of  these  estimates  is  discussed  below). 

There  are  two  reasons  we  consider  Eq.(6.51)  our  preferred  form  of  the  L-S  integral 
equation: 

(i)  As  shown  by  the  estimates  above,  only  two  nontrivial  terms  in  the  equation,  (6.51c) 
and  (6.51e),  remain  sizable  in  the  high  contrast  limit.  Both  of  these  terms  involve 
the  gradient  of  the  inverse  of  density,  and,  therefore,  give  rise  to  surface  contribution 
from  interfaces  of  regions  of  large  density  ratios. 
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(ii)  The  term  (6.51e)  involves  an  expression  proportional  to  the  stress  tensor  a, 


7h  6ij  diui  +  ( diuj  +  dj ut)  = 


An  °ij 


(6.52) 


Contraction  of  this  quantity  with  the  gradient  of  the  inverse  density,  i.e. ,  with  the 
normal  to  the  surface  of  density  discontinuity,  is  proportional  the  traction  vector, 
which  is  known  to  be  continuous  across  that  surface.  The  latter  property  facilitates 
discretization  of  the  equations  in  the  case  of  a  discontinuous  density  distribution, 
since  the  product  of  the  gradient  of  the  density  ratio  (a  surface  delta  function)  and  a 
continuous  function  is  unambiguously  defined. 


The  estimates  of  the  magnitudes  of  the  terms  in  Eq.(6.51)  follow  from  an  important 
scaling  property  of  Eqs.  (3.9c)  and  (6.51):  while  the  coefficients  involving  gradients  of  p^/p 
grow  proportionally  to  p/p^,  their  product  with  the  solution  remains  finite,  since  the  stress 
tensor  (6.52)  is  finite  in  the  high  density  limit,  and  thus  the  gradients  of  the  displacement 
in  Eq.(6.51)  are  small: 

4uj~^  (6.53) 

(we  recall  that  Eq.(3.10)  implies  ~  £  ~  1  for  finite  values  of  refraction  coefficients).  This 
is  the  reason  why,  in  Eq.(6.51),  the  terms  without  derivatives  of  the  density  and  factors  p/po, 
are  small.  In  particular,  the  fifth  term  (6.51e)  in  the  equation,  involving  the  gradient  of  £/( 
(and  thus  possibly  a  surface  delta  function  due  to  a  discontinuity  in  p)  is  small  compared 
to  the  dominant  terms,  (6.51c)  and  (6.51e). 


6.3.1  Matrix  elements  for  the  “basic”  form  of  second-order  equations:  general 
expressions 

The  most  basic  form  of  the  L-S  equation,  (6.48),  gives  rise  to  the  nontrivial  Galerkin  matrix 
elements1 

Aap  =  ~  f  d3ri  J  d3r2  C(r i)  {di  4  9(r i  -  r2))  (1  -  p{ r2))  Y^(r2))  ,  (6.54a) 

O  =  ~k~2  J  d3r1  f  d3r2C(ri)  (aJdUr,  -r2)) 

d2  [(!  -  A(r2))  4  V>J(r2)]  ,  (6.54b) 

AS  =  k~2  J  d3G  /  d3r2  V’a(r1)  {dldigir,  -r2)) 

4  [/i(r2)  {4  V$(r2)  +  4  ^(r2))]  (6.54c) 

1In  order  to  simplify  the  notation,  we  assume  here,  temporarily,  p0  =  \Q  =  1. 
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(the  other  matrix  elements  involve  single  integrals  only).  After  integrating  by  parts  with 
respect  to  r j .  we  find,  alternatively 


Aal  =  j  d3rl 

J  d3r2  (d\  COb))  (d[  g{ rx  -  r2))  (l  -  p( r2))  ^(r2))  , 

(6.55a) 

43  -  *-2  / 

d3r1  J  d3r2  (d\  (4  ~  r2)) 

4  [(!  -  A(r2))  ^^(r2)]  , 

(6.55b) 

4$  =  -  v2 

J  d3r1  J  d3r2  (9jC(r i))  (4#(ri  -  r2)) 

4  [Mr2)  ( 4  V’^(r2)  +  4  V^(r2))] 

(6.55c) 

CO 

r~0 

II 

S3 

J  d3r2  (d[  V’a(ri))  {d\  g{rl  -  r2))  (l  -  p(r2))  4?(r2))  , 

(6.56a) 

43  =  / 

d3r1  J  d3r2  (d[  ^fo))  (d\  g{r1  -  r2)) 

4  [(!  -  A(r2))  ^^(r2)]  , 

(6.56b) 

43  =  -  v2 

J  d3r:  J  d3r2  (4C(ri))  (4  4ri  -  r2)) 

4  [Mr2)  (4  V$(r2)  +  4  ^(r2))] 

(6.56c) 

Expected  features  of  the  matrix  elements.  We  note  that  in  Eqs.  (6.55)  the  basis 
function  appears  only  in  its  divergence.  For  linear  basis  functions  (Section  6.1.1),  this 
divergence  involves  constant  basis  functions  supported  on  tetrahedra  and  delta  functions 
supported  on  facets  located  on  the  boundary  dQ  of  the  object,  since  (as  discussed  following 
Eq.(6.3))  the  other  delta-function  contributions  cancel  pairwise,  due  to  continuity  of  the 
basis  function. 

Eqs.  (6.56)  involve,  instead  of  the  divergence  of  a  general  element  of  the  strain  tensor 
associated  with  this  basis  functions.  However,  the  regularity  structure  of  this  expression  is 
similar  to  the  previous  one. 

On  the  other  hand,  derivatives  of  the  basis  function  rtpa  do,  in  general,  include  delta- 
function  contributions  from  facets  shared  by  tetrahedra  supporting  the  basis  function  (since 
material  parameters  may  have  different  values  on  different  tetrahedra). 

Further,  Eqs.  (6.55b)  and  (6.55c),  as  well  as  Eqs.  (6.56b)  and  (6.56c),  include,  in  general, 
delta  functions  due  to  discontinuities  of  the  material  parameters  on  facets.  In  Eqs.  (6.56b) 
and  (6.56c)  the  product  of  these  delta  functions  and  the  derivatives  of  the  basis  function 
should  be  well  defined,  for  reasons  analogous  to  those  in  acoustics:  We  note  that  the  sum  of 
the  matrix  elements  (6.55b)  and  (6.55c)  the  considered  delta  functions  appear  only  through 
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the  expression 


J  d 3r2  (d[g(r  1  -r2)) 

[(4  A(r2))  di  'ipp(r2)  +  (4  g{r2))  (4  Y$(r2)  +  4  ^g(r2))] 

[  d2r2(d[g{r1-r2)) 

Jh 

nJ2  [(A(t2+)  -  A(*2-))  <4  4  ^(r2) 

+  (m(*2+)  -  M*2— ))  (4  V’^(r2)  +  4  ^(r2))] 


=  -  /  d2r2  (44ri  -r2)) 


n0 


'h 


{  [A(*2+)  <5y  4  ^(r2)  +  M*2+)  (4  ^J(r2)  +  4  ^(r2))] 
-  [(i2+  t2-)]  } 


=  -  /  d2r2  (4  s(ri  -  r2))  {er.- f(t2+)  -  a, .-(t2— )} 


'/2 


(6.57) 
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6.3.2  Matrix  elements  for  the  “high- contrast”  form  of  second-order  equations: 
general  expressions 

We  denote  in  the  following  by  the  subscripts  (a)  to  (f)  the  contributions  to  the  matrix 
elements  corresponding  to  terms  in  Eqs.  (6.51a)  to  (6.51e): 


S  =  J  d 3r 

V’aWV’M1’)  > 

(6.58a) 

(6)  -  A-2  / 

J 

d3r  (<9*  V4(r)) 

(!  -  CA(r))  <*ij  5fc^(r)  -  C^(r)  (5*  Y$(r)  +  4  Y^(r))  , 

(6.58b) 

(c)  -  A-2 

J 

^a(f)  5ij  5fc^(r)  +  ^(r)  (9V^(r)  +  , 

(6.58c) 

(d)  _  1.-2 

a/3  —  K 

J  d3rx  J  d3r2  (#1  ^^(rj)  5(ri  -  r2) 

4  (!  -  ^(r2))  sij  d|Y$(r2)  -  C/i(r2)  (4  V’^(r2)  +  4  ^(r2)) 

=  k~2  J  d3ri  J  d3r2  (^L(ri))  [d{g(r1  -  r2)) 

4  (!  -  ^(r2))  Sij  ^2^(r2)  -  ^(r2)  {4  V’Jfe)  +  4  >  (6.58d) 

A{;1  =  -  k-2  j d3r,  J d\  (a*  aj  *(: r,))  sK  -  r2)  («2^y) 

*7A(r2)  6ij  92^(r2)  +  ^(r2)  (d2^(r2)  +  d^(r2)) 

=  fe"2  J  d3rx  J  d3r2  (<9[  V’L(r1))  sK  -  r2))  (V2^_) 

?7A(r2)  Sij  92^(r2)  +  ^(r2)  (S^(r2)  +  •  (6-58e) 
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6.3.3  Matrix  elements  with  composite  linear  basis  functions 

After  carrying  out  the  index  algebra  in  Eqs.  (6.58),  the  matrix  elements  can  be  represents 
in  terms  of  the  scalar  basis  functions  as  follows: 


A —  A 

-iWc,m;v|8,n  —  umn 


d3r0VQ(r)</>V/3(r)  , 


(6.59a) 


4],m;v(iln  =  -k  2  /  d3r  (avva(r)) 


^Va,m;vfl,n  =  -  ^’2  /  d3r  0Va  (r)  (  9 


[?A(r)  5ni  5im  +  Mr)  (<5ni  Slm  +  Snm  8h)\  (dl  0V/3(r)) 
Po  \  P(r) 


A 


(d) 


va,m;v^,n 


=  k  2  I  d3r1  /  d3r2  ^  -  r2))  (  ^ 


P(r)/  Po 

[CA(r)  5nZ  5im  +  £/»  {sni  8lm  +  8nm  <JK)]  ( dl  0V/3( r))  , 

=  AT 2  J  d3^  j  d3r2  (^"^(rj))  -r2)) 

4  {  [(!  -  £A(r2))  ^nl  5ij  ~  ^(r2)  (<*m  <4  +  Snj  4*)]  ( 4  <t>w^2))  } 

=  -k^2  j  d3r1  J  d3r2  (3]71  <9(  0Va(ir))  (d{  g(r1  -  r2)) 

(!  -  ?A(r2))  Kl  Sij  -  ^(r2)  iSni  6lj  +  5nj  <*K)  ]  (4  0V/J(r2)) 

=  -k^2  f  d3rx  J  d3r2  (Sj77  ^(rq))  p(fl  -  r2) 

4  {  [(!  -  £A(r2))  ^nl  5ij  ~  Mr2)  4m  <4  +  Snj  ^ii)]  (4  0V/3(r2))  } 

Po 

!p(r2 


PA(r2)  ^nl  5ij  +  P4r2)  4ni  5lj  +  Snj  5li) ]  (4  ^sfo)) 


k  2  /  d3rx  J  d?r2{d^d\^a(rl))g{Yl 


-  r0 


a?_P0_ 

2p(r2 


PA(r2)  Snl  6ij  +  P^(r2)  4ni  6lj  +  <*k)  1  (4  </>V/3(r2)) 


(6.59b) 


(6.59c) 


(6.59d) 


(6.59e) 


In  the  following  we  specialize  to  the  case  of  material  parameters  constant  on  tetrahedra, 
and  to  piecewise  linear  basis  functions  0v(r).  The  matrix  elements  Eq.(6.59)  exhibit  then 
analogies  to  the  acoustic  problem  with  piecewise  linear  basis  functions,  since  a  derivative  of 
a  linear  basis  function  is  a  constant  -  plus,  possibly,  a  delta-function  at  the  boundary  of  the 
tetrahedron.  The  main  difficulty  lies  in  treating  the  possible  delta- function  contributions. 

Eqs.  (6.59d)  and  (6.59e)  appear  in  several  equivalent  forms  related  by  integration  by 
parts:  they  differ  in  the  numbers  of  derivatives  acting  on  the  Green  function  and  on  the 
basis  functions  ipa  or  (Eq.(6.9)).  The  forms  with  fewer  derivatives  acting  on  on  the 
Green  function  are  more  appropriate  for  computing  matrix  elements  at  small  distances 
(for  overlapping  basis  functions’  supports),  and  the  other  forms  are  more  useful  for  larger 
distances,  where  the  derivatives  of  the  Green  function  reflects  the  large-distance  behavior 
of  the  matrix  elements. 
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6.3.4  Matrix  elements  with  elementary  linear  basis  functions 

Eqs.  (6.59)  involve  “composite”  linear  scalar  basis  functions  associated  with  vertices,  i.e., 
supported  on  sets  of  tetrahedra.  In  the  following  we  use  the  relations  of  Section  6.1.1  to 
express  these  functions  in  terms  of  “elementary”  basis  functions  supported  on  individual 
tetrahedra  or  facets  (the  latter  arise  as  a  result  of  differentiation  of  the  volumetric  basis 
functions). 

We  give  below  explicit  expression  for  the  individual  contributions  to  the  matrix  elements. 
We  will  generally  concentrate  on  the  expressions  with  the  largest  numbers  of  derivatives  of 
the  basis  functions,  since  such  terms  may  contain  surface  contributions. 

The  term  A In  this  simplest  case  we  obtain  a  sum  of  contributions  of  tetrahedra  t 
shared  by  the  vertices  va  and  v^,  i.e.,  belonging  to  the  set  Ta  n  T p.  Each  term  is  an 
analytic  expression  for  the  integral  of  a  product  of  two  linear  functions  supported  on  the 
tetrahedron  t. 


The  terms  and  A^:  The  term  A^  involves  derivatives  of  the  basis  function  i/)Va,m 
which,  according  to  Eqs.  (6.1)  and  (6.3),  can  be  expressed  in  terms  of  basis  functions 
supported  on  tetrahedra  as 


s>v»  =  -  E 


I  t 

L  Vq, 


Ka,taXta{*)+  ^ 


«,/< 


tfn'fa.ta5! 


(r) 


fa&dtc 


(6.60) 


where,  we  recall,  hVa  t  the  height  of  the  tetrahedron  ta  measured  from  the  vertex  vQ:, 
hVy  t  is  the  exterior  unit  normal  to  the  tetrahedron  facet  opposite  the  same  vertex,  and 
ta  is  the  unit  normal  to  the  face  fa,  in  the  direction  exterior  to  the  tetrahedron  ta. 

The  outer  sum  in  Eq.(6.60)  runs  over  tetrahedra  ta  £  Ta,  sharing  the  selected  vertex 
va.  The  inner  sum  runs  over  faces  fa  £  dta  of  each  tetrahedron  ta.  As  already  discussed 
in  Section  6.1.1,  in  the  latter  sum  the  delta- function  contributions  vanish  on  the  outer 
boundary  of  the  set  dT a  and  cancel  pairwise  for  interior  facets  fa,  although  not  for  the 
boundary  facets  (/Q  £  dfl).  Therefore,  Eq.(6.60)  reduces  to 


d>vQ(r)  = 


-  T  — 

hv  t 

f  (Z'T  Vck,  t-Q 
lOL  C  J.  ck 


n. 


Xt„(r) 


-  E 

fa&FaPdCl 


Va,/, 


.CO 


(6.61) 


where  T a  is  the  set  of  facets  sharing  the  vertex  vQ. 

The  term  A ^  involves  only  the  basis  function  i/)Va,m,  he.,  according  to  Eqs.  (6.1)  and 
(6.2),  is  given  by 


0v„(r) 


n 


lva,t0 


Va,ta 


■  (r 


Xta(r)  , 


(6.62) 


where  (r)  is  the  characteristic  function  of  the  tetrahedron  ta.  This  basis  function  is 
multiplied  by  the  gradient  of  the  inverse  density,  which  is  proportional  to  a  surface  delta 
function  on  each  facet  on  which  the  density  is  discontinuous.  In  the  matrix  element  A ^ 
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(6.63) 


such  delta- function  contributions  arise  on  facets  shared  by  the  tetrahedra  ta  and  tp  (the 
tetrahedra  may  be  identical). 

In  both  of  the  terms  A ^  and  A ^  the  basis  function  <j>  appears  only  through  the 
expression 

Sim(^vJ8,n(r))  :=  (!  -  £A(r))  5im  ^  “  ^(r)  ^vp,n(r)) 

=  [(!  -  £a(F))  Kl  5im  ~  £/»  iSni  6lm  +  5nm  <^)  ]  Kp(r)) 

(see  Eq.(6.10)),  which  is  related  to  the  stress  tensor  atm  associated  with  the  basis  function 

^V/3,n(r)- 

Eq.(6.63)  has  to  be  represented,  by  using  Eqs.  (6.1),  (6.2),  and  (6.3),  in  terms  of  the 
basis  functions  cj)v  t  supported  on  the  individual  tetrahedra.  The  result  is 

Stm(^,n(  r))=  E  Ximnl(tp)dl(l)Vl3,tl3(r) 

tp^Tp 


nv 


h  '  ~  V/3 

L  lvp,tp 


,tpXt0(r)+  £ 


fp&dt.p 


(6.64) 


where 


Ximnl(^(s)  ■—  (1  Ca(^))  ^ im  ^nZ  +  £/*(*/?)  (^m  ^mZ  +  ^ il  ^ ran )  r  (6.65) 

and  £\(tp)  and  (t;3)  are  the  values  of  the  coefficients  and  on  the  tetrahedron  tp.  It 
is  worth  noting  here  that,  by  construction,  the  coefficients  Ximni  have  the  same  symmetry 
properties  (Eq.(3.12))  as  the  general  elasticity  tensor  for  an  anisotropic  medium 


X  =  X 

imnl  mini  ’ 


imnl 


=  x 


imln  ’ 


v  _  v 
yv  imnl  nlim 


(6.66) 


We  will  use  this  property  in  the  following.  A  related  remark  is  that  the  manipulation  we 
are  carrying  out  now  could  be  probably  generalized  to  anisotropic  media,  and  the  resulting 
stiffness  matrix  would  have  a  similar  structure  to  that  in  the  present  case. 

By  changing  the  order  of  summation,  Eq.(6.64)  can  be  recast  in  the  form 


‘='im(^,v/3,ra(r)) 

=  ~  E  Ximnl(tp)f^—iilv0,tpXtp(r) 

-  E  K,fp(r)6fp(r)  E  XimnAtp)™lU,tp  , 

fl 3£Ts  t0er//3 


(6.67) 


where  T ^  is  the  set  of  (up  to  two)  tetrahedra  sharing  the  facet  f  g.  The  last  sum  in  Eq.(6.67) 
is  thus,  in  general,  a  difference  of  two  terms  associated  with  the  tetrahedra  adjacent  to  the 
given  face  fg.  Therefore,  we  can  write  Eq.(6.67)  as 


1—1  im  vp  ,n  (t)  ) 

=  "  E  Ximnl(^)u - hvp,tpXtp( r) 

tp&Tp  v/3’*/3 

-  E  [Ximnl(tf3+)  ~  Xim.nl(t0-)]  ^fp  6 fp(r)  <f>vp,fp(r) 
fpzFpndn 


(6.68) 
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where  tg+  and  t,g—  are  tetrahedra  on  the  positive  and  negative  sides  of  the  facet  fg  (relative 
to  the  direction  of  the  normal  ny  to  the  face).  We  also  used  here  the  fact  (discussed  in 
connection  with  Eq.(6.61))  that  the  delta-function  contributions  cancel  pairwise  on  all  facets 
fp  except  those  located  on  the  exterior  object  boundary  <917. 

Finally,  we  note  that  we  can  drop  the  delta-function  contributions  (i.e. ,  the  second  sum 
in  Eq.(6.69)),  altogether:  the  reason  is  that  we  can  think  of  the  object  boundary  <917  as  an 
interface  between  the  object  and  the  external  region  discretized  with  “fictitious”  tetrahedra 
having  properties  of  the  background  medium  (£A  =  1  and  =  0,  hence  X  =  0).  With  this 
interpretation,  the  surface  <917  is  no  longer  a  boundary  (we  have  used  a  similar  argument  in 
acoustics).  Hence,  we  have,  eventually, 

Him(^,n(r))  =~  J2  r)  ■  (6-69) 


It  will  be  also  useful  to  introduce  a  quantity  related  directly  to  the  stress  tensor, 


Sim(^vj8,n(r))  :=  , 

which  can  be  represented,  in  analogy  to  Eq.(6.69),  as 
S*m(V>V/3>n(r)) 

=  “  2  r) 
tpeTp 

-  [Elmnlit0+)  ~  Eimnl(t0-)\  4/3  4/3  (r) 


(6.70) 


(6.71) 


where 


Eimni(t0)  :=  Vx(tg)  6im  5nl  +  rj g)  ( 8in  5ml  +  Su  5mn)  ,  (6.72) 


and 


V\(tg)  ■= 


P(tp) 

Po 
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pits) ,  u  ,  p(tp) 

Po  P  Ao 


(6.73) 


(cf.  Eqs.  (3.10)  and  (3.11)).  The  expression  Him(rfV0^n(r))  appears  in  the  matrix  element 

A^. 

In  terms  of  the  tensors  (6.63)  and  (6.70),  the  contributions  to  the  matrix  elements  are 
then  given  simply  by 


and 


Avltm;vp,n  =  -k  2  J  d3r  (<9>Va(r))  Him(^Vj8)n(r))  (6.74) 

A(flm]V0,n  =  -  k2  J  d3r  (j)Vct  (r)  ■  (6-75) 
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The  terms  and  A^:  The  remaining  terms  in  the  matrix  elements  involve  several 
common  expressions: 

Second  derivatives  of  the  basis  function  (j)a.  These  expressions  (Eq.(6.8))  were 
given  in  Section  6.1.1.  Clearly,  they  cannot  be  used  directly  in  evaluating  the  matrix 
elements:  derivatives  of  the  delta  function  have  to  be  eliminated  in  favor  of  derivatives  of 
the  Green  function.  We  will  carry  out  this  rearrangement  in  the  following. 


The  stress  tensor  associated  with  the  basis  function  (j)V0.  In  Eqs.  (6.59d)  and 
(6.59e)  the  basis  function  d>Vfl  appears  only  through  the  quantity 


3v(^V/S,n(r))  :=  £A(r)  6ij  &  ^ ,n(r)  +  C/i(r)  (5*  lKf,,n(r)  +  &  ^.nW) 

"(1  -  CA(r))  ^  -  e„(r)  (sni  5tj  +  *BJ.  <y  1  ( dl  ^(r)) 


(6.76) 


(cf.  Eq.(6.63)),  related  to  the  stress  tensor  aV}  associated  with  the  basis  function  r/>v^  n(r). 
It  can  be  expressed,  as  in  Eq.(6.69),  in  terms  of  constant  basis  functions  supported  on  the 
tetrahedra. 

Similarly  to  Eqs.  (6.74)  and  (6.75),  Eqs.  (6.59d)  and  (6.59e)  take  the  form 
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(6.78) 


Full  expressions  for  and  A^e\  After  substituting  Eq.(6.8)  in  Eq.(6.77)  we 

integrate  one  of  the  terms  by  parts  in  order  to  eliminate  the  derivatives  of  the  delta  function 
in  the  expression  for  the  derivatives  of  the  basis  function  <f)Va-  We  obtain,  instead,  an 
additional  derivative  of  the  Green  function  contracted  with  the  normals  to  the  facets  on 
which  the  derivative  of  <pV;y  is  supported  (hja  c ^  g(r1  —  r2)  in  the  third  term  of  Eq.(6.79) 
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below).  The  result  is 
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(6.79) 


Similarly,  we  obtain 
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Finally,  by  substituting  Eq.(6.69)  in  Eq.(6.79),  we  obtain  a  sum  of  several  terms,  involv- 


47 


ing  integrals  over  tetrahedra  and  facets, 
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Similarly,  Eq.(6.71)  substituted  in  Eq.(6.80)  yields 
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(6.82) 
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By  reversing  the  order  of  summation,  the  above  expression  can  be  written  as 
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All  the  contributions  to  the  r2-integrals  in  Eq.(6.82)  come  from  discontinuities  of  the 
inverse  density,  Po/p(//j+)  —  Po/p(f/3~ )>  where  p(//3+)  and  p(fp-)  are  densities  in  tetra- 
hedra  on  the  positive  and  negative  sides  of  the  face  fp,  as  defined  by  the  normal  t 
(which  is,  by  definition,  exterior  to  the  tetrahedron  tp). 

In  the  alternative  expression  Eq.(6.83)  come  from  discontinuities  of  the  inverse  density, 


6.3.5  Summary  of  the  expressions  for  the  “basic”  matrix  elements 

The  fundamental  matrix  elements,  required  in  the  computation  of  the  matrix  elements 
the  second-order  formulation  are  listed  in  the  following.  Symbols  (d)  and  (e)  indicate  the 
contribution  in  which  the  given  matrix  element  appears. 


tetrahedron-facet  matrix  elements: 
constant-constant: 

TF2(d)  Al{ta-Jp)=f  d3r1  f  d2r2  dl  g(r1  -  r2) 
Jta  J  f  p 

TF3(d)  Aln{taJ p)  =  h3fp  f  d3r1  f  d2r2  dl  &>  g{r1 

J  toL  J  f(3 


-  rn 
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facet-facet  matrix  elements: 


constant -  const  ant : 

FFl(e) 

A(fa‘Jp)=  f  d2 rl 
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The  subscript  “n”  in  (•)  is  used  to  indicate  matrix  elements  involving  the  normal  deriva¬ 
tive  of  the  Green  function. 


6.4  Representation  of  matrix  elements 

We  recast  here  matrix  elements  evaluated  in  the  previous  Subsections  in  a  general  represen¬ 
tation  used  in  the  following  in  computation  of  the  compressed  stiffness  matrix.  This  purpose 
of  using  this  representation  is  to  facilitate  programming  and  allow  an  efficient  ordering  of 
operations  in  the  matrix  construction  stage  (as  described  in  Section  6.6). 

We  give  here  a  somewhat  schematic  expression  for  the  matrix  -  the  details  of  its  imple¬ 
mentation  will  be  given  in  Section  6.6. 

The  representation  has  an  almost  factorized  form2 
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In  this  expression: 

1.  The  indices  a  and  (3  are  unknown  numbers. 

2.  The  indices  Vj,  j  =  1,2,  label  specific  contributions  to  the  matrix.  In  the  practical 
implementation  they  are  simply  integers  uniquely  identifying  the  given  contribution 
to  the  matrix  and  the  given  coefficient  C. 

3.  The  matrix  B  is  a  set  of  coefficients  0  or  1  indicating  which  contributions  are  included 
in  the  sum. 

4.  The  multi-indices  Mj,  j  =  1,2  (generally,  M )  label  the  “fundamental”  basis  functions 
cf>M.  In  the  present  case,  they  are  constant  or  linear  basis  functions  supported  on 
tetrahedra  or  facets. 

5.  The  operators  V11 12  are  differential  operators  representing  partial  derivatives  specified 
by  the  (multi-)indices  /,  acting  on  the  scalar  Green  function  g.  In  practice,  derivatives 
only  up  to  the  second  order  will  appear. 

2Full  factorization  does  not  seem  possible,  because  some  matrix  elements  are  modified  by  hand  through 
integration  by  parts  and  similar  manipulations. 
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6.  The  coefficients  C^Vj^Ij  provide  a  link  between  the  unknowns  and  the  basis  functions. 
The  indices  Vj  are  additional  labels  indexing  possible  types  of  unknowns.  In  general, 
the  coefficients  C  (especially  C will  depend  on  the  material  parameters.  In  the 
following  we  consider  a  given  set  of  coefficients  as  a  sparse  matrix  whose  rows 

are  labeled  by  the  unknowns  a. 

In  addition  to  the  coefficients  C  (which  provide  a  mapping  from  the  unknowns  to  the 
fundamental  basis  functions),  we  also  need  an  “inverse”  mapping  D  from  the  fundamental 
basis  functions  to  the  unknowns.  As  we  discuss  later,  these  mapping  are  used  in  the  matrix 
assembly. 

The  idea  of  the  representation  (6.84)  is  to  provide  a  universal  scheme  for  matrix  element 
computation,  and  allow  full  flexibility  in  defining  and  modifying  the  form  of  the  integral 
equations. 

Sets  of  the  coefficients  C  are  precomputed  and  then  supplied  as  arguments  to  an  “assem¬ 
bly”  routine,  which  constructs  contributions  to  the  output  matrix  elements  from  the  “basic” 
matrix  elements  defined  in  terms  of  the  “fundamental”  basis  functions  and  the  Green  func¬ 
tion  or  its  derivatives.  The  same  coefficients  are  used  in  evaluation  of  the  MoM-Cartesian 
mapping  coefficients  V,  as  the  mapping  of  unknowns  will  be  expressed  in  terms  of  mappings 
of  the  fundamental  basis  functions. 

The  coefficients  C  are  also  incorporates  as  a  component  in  the  compressed  stiffness 
matrix.  The  reason  is  that  elements  of  the  map  VQ  for  the  unknowns  may  be  linear  com¬ 
binations  of  many  map  elements  for  the  fundamental  basis  functions;  therefore,  it  is  more 
economical  to  store  those  basic  mappings  and  convert  the  MoM  variables  into  equivalent 
sources  (and  vice  versa)  in  two  stages. 

In  the  above  description  we  assumed  that  the  far-held  part  of  the  matrix  is  represented 
not  only  in  terms  of  the  Green  function,  but  also  of  its  derivatives.  However,  this  type  of 
representation  has  to  be  carefully  examined  on  the  case-by-case  basis,  particularly  from  the 
point  of  view  of  storage.  If  we  decide  to  use  the  Green  function  only,  the  far-held  coefficients 
will  be  different  than  the  coefficients  C  above;  nevertheless,  storage  of  the  expansion  of  the 
mappings  for  the  fundamental  basis  functions  may  be  still  more  economical  than  for  the 
unknowns. 

In  the  expressions  below  the  symbols  C (t)  and  C (/)  refer  to  functions  constant  on 
tetrahedra  and  facets,  and  L(v,  t)  and  L(v,  /)  to  functions  linear  on  tetrahedra  and  facets; 
in  the  latter  case,  v  is  the  selected  vertex  (at  which  the  basis  function  has  value  1).  Extra 
overall  coefficients  are  included  in  the  B. 

6.4.1  Matrix  elements  for  second-order  equations 

We  give  here  examples  of  expressions  for  the  matrix  elements  in  the  second-order  formulation 
of  Eqs.  (6.74),  (6.75),  (6.81),  (6.82)  (these  are  the  most  complicated  ones). 

We  give  alternative  expressions  for  some  coefficients  C,  differing  by  the  presence  of 
absence  of  one  factor  n  ;  those  with  the  missing  factor  are  used  with  basic  matrix  elements 
involving  the  normal  derivative  of  the  Green  function.  For  notational  simplicity,  some  of 
the  constant  coefficients  are  not  displayed  in  these  formulae,  and  are  absorbed  in  the  factors 
B  appearing  in  Eq.(6.84). 
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The  contribution  A 
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In  terms  of  these  coefficients,  the  matrix  element  (6.81), 


fa  e  -^o  7 

(6.85a) 

fa  €  ^Ta  , 

(6.85b) 

fa  g  dTa  n  30  , 

(6.85c) 

fa  e  dTa  n  > 

(6.85d) 

t/3  G  ^  > 

(6.85e) 

(6.85f) 

written  in  the  form  of  Eq.(6.84) , 
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(6.86) 
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The  contribution  A (e); 
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In  terms  of  these  coefficients,  the  matrix  element  (6.82),  written  in  the  form  of  Eq.(6.84), 
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(6.88) 


6.5  Tetrahedron-tetrahedron  contributions  to  stiffness  matrix  elements 

Each  of  the  matrix  elements  for  “composite”  basis  functions  is,  generally,  a  sum  of  many 
contributions  of  pairs  of  “fundamental”  basis  function.  This  is,  in  particular,  the  case  for 
node  based  functions,  where  a  vertex  may  be  shared  by  many  tetrahedra  and  facets. 
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We  now  estimate  the  number  of  contributing  tetrahedra  and  facets,  e.g.,  the  number 
of  ta,  fa,  etc.,  entries  in  each  “row”  of  the  coefficient  matrices  of  Eqs.  (6.85)  and  (6.87). 
We  are  considering  here  the  generic  case  of  a  vertex  located  in  the  interior  of  a  regular 
tetrahedral  mesh  (approximately  constant  edge  lengths). 

First,  we  note  that  the  number  of  tetrahedra  ta  in  the  set  Ta,  or,  equivalently,  the 
number  of  facets  fa  in  the  set  dTa,  is  approximately  equal  to  the  number  of  triangles 
in  a  regular  triangulation  of  the  unit  sphere.  Since  the  area  of  a  triangle  is  A  ~  \/3/4, 
the  number  of  such  triangles  covering  the  sphere  is  about  4ir/A  ~  47r/(v/3/4)  =  167t/\/3. 
Secondly,  the  number  of  facets  fa  sharing  the  given  vertex,  i.e.,  the  number  of  elements  in 
the  set  Jra,  is  equal  to  the  number  of  edges  in  the  considered  triangulation  of  the  sphere, 
which  is  about  3/2  the  number  of  facets,  hence  (3/2)  167t/\/3  =  8\/3tt-  Hence, 

#Ta  =  #8Ta  -  -  29.0  (6.89a) 

and 

~  8\/3vr  ~  43.5  .  (6.89b) 

6.6  Construction  of  the  stiffness  matrix 

In  designing  the  solver  code,  we  have  to  allow  sufficient  flexibility  in  order  to  accommodate 
various  types  of  equations,  various  sets  of  unknowns,  and  expressions  for  matrix  elements. 


General  structure  of  matrix  elements.  All  the  expressions  for  the  matrix  elements 
involving  the  Green  function,3  i.e.,  Eqs.  (6.81),  (6.82),  and  (6.39)  -  (6.43)  have  the  struc¬ 
ture  indicated  in  Eq.(6.84).  More  explicitly,  for  a  given  equation  type,  a  matrix  block 
corresponding  to  two  fixed  types  of  unknowns4  has  the  form 


4*  =  E  E  E  E  ■ 

vi  V1  91,  bi  <?2,02  l 

(6.90) 

The  notation  here  is  as  follows: 


1.  Within  the  considered  block,  unknowns  are  labeled  locally  by  integer  indices  a,  (3. 

2.  The  sets  of  coefficients  C  are  stored  as  arrays  indexed  by  the  indices  Up  j  =  1,2;  the 
order  of  the  entries  in  the  array  is  arbitrary. 

3.  Similarly,  the  “basic”  matrix  elements  \-)  are  labeled  by  an  index  v,  which  defines 
the  type  of  the  matrix  element  (e.g.,  between  a  linear  function  on  a  tetrahedron  and 
a  constant  function  on  a  facet,  and  with  the  normal  derivative  of  the  Green  function). 
Actually,  in  the  actual  code  implementation  the  matrix  elements  A  are  not  stored, 
but  rather  computed  as  needed  (in  a  loop  which  ensures  there  no  replication  of  work); 
in  this  case  the  index  v  is  an  argument  of  the  routine  computing  the  matrix  elements 
and  specifies  which  element  is  to  be  evaluated. 

3The  matrix  elements  given  by  single  integrals,  as  well  as  projections  of  the  incident  or  scattered  waves 
on  the  basis  function,  can  be  treated  as  special  cases. 

4E.g.,  the  row  indices  a  may  refer  to  displacement,  and  the  column  indices  f3  to  pressure. 
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4.  The  indices  i/j,  j  =  1,2,  refer  to  consecutive  sets  of  the  coefficients  C;  they  form  an 
array  of  coefficient  sets,  and  may  be  ordered  in  an  arbitrary  way. 

5.  Accordingly,  the  coefficients  B(vl,v2)  define  weights  with  which  pairs  of  coefficients 
C  contribute  to  the  sum. 

6.  Similarly,  the  integer  indices  M(v1,v2)  specify  which  type  of  a  matrix  element  is  to 
be  used  with  the  given  coefficients  C . 

7.  The  indices  g1  and  g2  refer  to  g-element  types  (in  the  present  problems,  tetrahedra 
or  triangles),  and  the  indices  b1  and  b2  to  consecutive  numbers  of  basis  functions 
associated  with  a  given  g-element;  e.g.,  if  the  considered  basis  functions  are  linear  and 
the  g-element  is  a  tetrahedron,  the  index  b  may  take  values  from  1  to  4,  and  indicate 
the  tetrahedron  vertex  at  which  the  basis  function  value  is  1. 

8.  The  coefficients  provide  relations  between  the  unknowns  and  the  basis  func- 

a,  gj ,  Oj 

tions;  they  depend,  in  general,  on  the  material  properties  of  the  object. 

9.  The  quantities  A^I(g1,  b1;g2,  b2)  are  the  “basic”  matrix  elements,  such  as  listed  in 
Section  6.3.5.  As  we  mentioned  before,  the  index  v  may  be  used  to  specify  various 
types  of  matrix  elements. 

10.  The  integer  indices  Ij  in  the  coefficients  specify  vector  or  tensor  indices.  If 

C  is  a  rank-1  tensor  and  is  labeled  by  a  single  index,  I  j  may  be  identified  with  that 
index.  If  C  is  a  rank-2  tensor,  then  the  index  /  ■  may  refer  to  pairs  of  Cartesian 
indices;  for  example,  if  the  tensor  is  symmetric,5  the  consecutive  values  1,2,  ...  ,6  of 
Ij  may  correspond  to  pairs  (1, 1),  (1, 2),  ...  ,  (3,  3). 

11.  The  indices  I  in  the  matrix  elements  have  the  same  meaning  as  in  the  coef¬ 

ficients  C.  (In  practice,  I  will  be  usually  just  a  single  Cartesian  index.) 

12.  The  tensor  algebra  (index  contraction)  is  handled  by  taking  the  sum  over  £:  the  arrays 
I j{£)  and  I(t)  define  then  the  matching  tensor  indices,  their  pairs,  etc. 

Matrix  storage.  In  practice,  for  larger  problems  it  is  always  necessary  to  use  matrix 
compression.  Hence,  the  matrix  blocks  AU1 V2  are  stored  as  sparse  matrices,  involving  only 
couplings  between  spatially  close  geometry  elements.  Therefore,  an  important  aspect  of 
matrix  construction  is  determination  of  its  sparsity  structure,  based  on  distances  between 
the  geometry  elements;  the  problem  is  relatively  complex,  since  different  near-field  ranges 
will  be  needed  for  various  types  of  matrix  elements  (near-field  matrix  element  computation 
and  storage  is  expensive,  hence  we  want  to  reduce  the  near  field  range  as  much  as  possible). 
In  addition,  of  course,  we  will  have  to  store  additional  matrix  data  required  in  compression, 
such  as  coefficients  for  mapping  basis  functions  into  equivalent  Cartesian-grid  sources. 

Further,  parallel  implementation  of  the  solver  requires  utilization  of  geometry  segments 
(slices  or  stacks  of  slices)  rather  than  a  single  global  geometry.  Thus,  matrix  blocks  are 
always  constructed,  on  a  given  processor,  by  using  locally  available  geometry  segments. 

5The  symmetry  may  follow,  e.g.,  from  the  relations  (6.66). 
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Organizing  matrix  construction.  Obviously,  the  most  time  consuming  element  in  the 
matrix  computation  is  the  evaluation  of  the  basic  matrix  elements  A(g1,  b1;g2,  b2)  them¬ 
selves.  Therefore,  in  order  to  minimize  the  computational  cost  without  introducing  addi¬ 
tional  storage,  many  (perhaps  all)  matrix  blocks  should  be  created  incrementally,  but  in 
parallel  -  in  the  sense  that,  having  computed  a  given  matrix  element  A(g1,  b1;g2 ,  b2),  one 
should  add  its  contribution  to  all  relevant  matrix  elements.  One  sould  also  carry  out  far- 
held  subtractions  on-the-hy,  while  storing  near-held  matrix  elements  (this  is  particularly 
important  in  the  parallel  computation). 

The  envisaged  scheme  is  as  follows 

1.  Compute  blocks  of  coefficients  C^Uj\ 

2.  Compute  MoM-Cartesian  mappings. 

3.  Compute  maps  D  “inverse”  to  the  maps  C,  in  the  sense  of  mapping  a  given  geometry 
element  onto  a  set  of  unknowns.  We  dehne  separate  maps  D  for  the  types  of  geometry 
elements  involved  (here,  tetrahedra  and  facets)  and  for  the  types  of  unknowns  (e.g., 
pressure  and  displacement).  I.e.,  for  each  type  of  the  g-element  g  and  each  type  of 
unknown  u  we  dehne  a  matrix  D^-a\  in  integer  compressed  sparse  row  format  (integer 
CSR);  rows  of  this  matrix  will  be  labeled  by  the  g-elements  and  each  row  will  contain 
a  list  of  unknown  numbers  corresponding  to  the  basis  function  numbers  b. 

4.  Compute  sparsity  patterns  of  matrix  blocks  (or  “sub-blocks”,  if  various  types  of  ele¬ 
ments  are  involved  -  e.g.,  interface  and  interior  elements,  etc.).  In  the  parallel  code, 
the  computed  sparsity  patterns  must  take  into  account  location  of  elements  in  geom¬ 
etry  segments;  e.g.,  row  indices  must  always  refer  to  the  “home”  geometry  slice,  while 
column  indices  may  refer  to  elements  in  the  stack  of  slices. 

5.  On  the  basis  of  nonzero  coefficients  and  the  sparsity  patterns,  create  lists  of 

geometry  elements  that  will  be  required  in  the  matrix  computation  (or,  maybe,  set 
hags  for  elements  in  the  input  geometries).  For  instance,  we  may  use  all  tetrahedra 
in  the  geometry,  but  only  some  triangles  (those  at  interfaces),  etc. 

6.  Carry  out  loops  over  selected  geometry  elements.  For  each  pair  of  elements: 

(a)  compute  all  basic  matrix  elements  (as  required  by  the  coefficients  and 

sparsity  patterns); 

(b)  compute  far-held  subtractions  (in  terms  of  the  MoM-Cartesian  mappings)  and 
modify  the  evaluated  matrix  elements; 

(c)  add  their  contributions  to  the  required  blocks. 

Schematically,  if  we  assume  the  sparsity  patterns  are  stored  in  terms  of  g- 
elements,  the  structure  of  the  loop  is 
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for  ngl  =  1,  . . .  ,  nngl 

/ /  get  geometry  data  7j  of  g-element  ngl 

//  compute  and  store  MoM-Cartesian  mapping  coefficients  for  g-elem.  ngl: 

W1  =  W{1\  71) 

//  get  number  of  near  elements  in  the  sparsity  pattern: 

nng2  =  . . . 

for  ng2  =  1,  . . .  ,  nng2 

//  get  geometry  data  72  of  g-element  ng2 

/ /  compute  MoM-Cartesian  mapping  coefficients  for  g-elem.  ng2: 

W2  =  W^(  72) 

//  store  W 2  only  for  j2  in  home  slice 

//  get  basic  MoM  matrix  element  for  the  g-elements  j2: 

AM°M  =  A(  71,72) 

//  get  basic  far- field  matrix  element  for  the  g-elements  7l5  72: 

/ /  [this  should  be  modified  to  take  into  account  derivatives  of  g] 

AFal  =  Wl9Wj 

//  get  numbers  of  unknowns  associated  with  71  and  72: 

kkxl  =  A'(D(1)(ngl))  (6.91) 

kkx2  =  A'(D*'2')(ng2)) 
for  kxl  =  1,  . . .  ,  kkxl 

//  get  unknown  number  from  inverse  mapping: 
nxl  =  (ngl,  kxl) 

//  get  the  coefficient  relating  the  unkn.  nxl  and  the  g-elem.  ngl: 
cx  =  C ^  (nxl,  ngl) 
for  kx2  =  1,  . . .  ,  kkx2 

//  get  unknown  number  from  inverse  mapping: 
nx2  =  Z?^(ng2,  kx2) 

//  get  the  coefficient  relating  the  unkn.  nx2  and  the  g-elem.  ng2: 
c2  =  C'^2^(nx2,ng2) 

//  add  computed  contribution  to  the  output  matrix  element: 

A(nxl,nx2)  +=  <7  c2  (AMoM  -  ^lFar) 
endfor  kx2 
endfor  kxl 
endfor  ng2 
endfor  ngl 

We  stress  that  in  the  above  scheme  the  MoM-Cartesian  mapping  is  computed  on-the- 
fly.  Out  of  all  the  mapping  coefficient  sets  XV  we  store  only  those  which  are  needed  in  the 
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matrix- vector  multiplication . 


7  Implementation  of  the  volumentric  integral-equations  code 
for  elasticity 

We  describe  here  the  initiated  implementation  of  the  integral  solver  code  for  elasticity.  The 
work  is  in  progress,  and  its  main  goal  is  to  develop  a  comprehensive  library  of  functions 
covering  various  types  of  integral  equations  described  in  this  report. 

7.1  The  code  structure 

The  present  version  of  the  code  contains  several  routines  written  to  allow  testing  of  various 
types  of  integral  equations  in  acoustics. 

The  general  idea  of  the  code  is  to  use  combinations  of  as  simple  routines  as  possible. 
In  particular,  geometry  data,  mapping  between  geometry  elements  and  unknowns,  and 
material  data  are  kept  separate,  and  handled,  as  far  as  possible,  by  separate  routines. 
Thus,  for  example,  there  are  separate  routines  for 

1.  Computing  mapping  between  geometry  elements  ( “g-elenrents” )  and  unknowns,  using 
as  input  geometry  and  material  data. 

2.  Generating  sparsity  pattern  of  the  near- field  matrix  blocks,  using  as  input  geometry, 
geometry-unknown  mapping,  and  compression  data. 

3.  Generating  sparsity  pattern  of  the  near- field  matrix  blocks,  using  as  input  geometry, 
geometry-unknown  mapping,  and  compression  data. 

4.  Generating  blocks  of  “basic”  matrix  elements  between  the  “fundamental”  basis  func¬ 
tions  supported  on  g-elenrents,  using  as  input  geometry,  geometry-unknown  mapping, 
and  matrix  sparsity  patterns. 

5.  Generating  sets  of  material-depending  coefficients  appearing  in  matrix  elements  for 
the  actual  equations  (defined  in  terms  of  the  “composite”  basis  functions). 

6.  Assembling  matrix  blocks  of  the  final  matrix  for  the  given  integral  equation,  using 
as  input  geometry  (its  connectivity  data),  matrix  blocks  for  the  “fundamental”  basis 
functions,  and  the  material-dependent  coefficients  computed  in  the  previous  step. 

7.  Generating  MoM-to-Cartesian  mapping  coefficients  for  the  “fundamental”  basis  func¬ 
tions,  using  as  input  geometry  and  geometry-unknown  mapping. 

8.  Assembling  MoM-to-Cartesian  mapping  coefficients  for  the  “composite”  basis  func¬ 
tions,  using  as  input  geometry,  the  previously  computed  mapping,  and  the  material- 
dependent  coefficients. 
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Blocks  of  the  MoM  matrix  are  computed  according  to  the  formula 


A 


(") 

a1=K1(i/,K.)  p1(iy,K,r1)-j-I1(iy,K,t ) ,  <y.2=K2(y,n)  p2(v,K,,r2)+I2(i',K,t) 

^max  K-u 

+  =  EESEEE»t)ESE 

a i  C2  71  72  i/=l  k=1  r\  t=  1 

C'(1)(^, «,  ruA(zy>  *))  «,  r2,I2(u,  re,  t))  , 


(7.1) 


which  is  more  explicit  from  of  Eq.(6.90). 

The  expression  (7.1)  represents  a  set  of  ^max  matrix  blocks,  z/  =  1,  ...  ,  z/max .  The  sum 
on  its  r.h.s.  has  to  be  understood  as  a  loop  in  which  contributions  to  the  individual  matrix 
blocks  are  being  accumulated.  In  this  loop  the  unknown  numbers  a1  and  a2  are  expressed 
in  terms  of  the  tensor  indices  /1  and  I2,  which  are  given  by  the  tensor  multiplication  table 
for  the  given  matrix  block  number  v  and  the  given  pair  k  of  mapping  data.  The  tensor 
index  I  -  selecting  an  element  from  the  set  of  the  basic  matrix  elements  -  is  also  determined 
by  the  tensor  multiplication  table.  The  unknowns  indices  oq  and  a2,  labeling  the  output 
matrix  element,  are  functions  of  the  sum  (loop)  indices  za  k.  r1;  r2,  and  t. 

The  actual  implementation  of  the  expression  (7.1),  in  terms  of  loops  and  variables  in 
the  code,  is  given  by  Eq.(7.18).  To  facilitate  the  comparison  of  the  formula  (7.1)  and  the 
pseudo-code,  we  show  here  the  correspondence  between  the  relevant  quantities,  particularly 
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the  indices  (labeled  by  j  =  1,2): 


^max  “»•  nnu 

o-j  -»•  kgt  j 
lj  ->•  ngj 
z/  — ►  mi 
AT,,  — ►  kkcc 

=  ppcb  [nu] ->pcob->na 
k  — ►  kcc 

^  -  cb 

=  ppcb  [nu] ->pcob->paa [kcc] 

— >  ncaj 

=  ppcb  [nu] ->pcob->pia [kcc] 
c2  — >  ncaj 

=  ppcb  [nu] ->pcob->pja[kcc] 

K  ■  — >  kkxcj 

=  ppcaj [ncaj]  ->kkx 
r  •  — >  narj 
Pj  ->•  nr  j 

=  ppcaj  [ncaj]  ->pcsr->p j a [narj] 
t  — >  nat 
/j  — ►  nil 

=  ppcb  [mi] ->ppcot [kkc] ->pia [nat] 
I2  —■ >  ni2 

=  ppcb  [nu] ->ppcot [kkc] ->pja [nat] 
/  — >  nia 

=  ppcb  [nu] ->ppcot [kkc] ->paa [nat] 

a  ->  nxj 


(number  of  matrix  blocks)  , 

(g-element  type)  , 

(g-element  number)  , 

(matrix  block  number)  , 

(number  of  pairs  of  mapping  data  sets)  , 

(index  of  a  pair  of  mapping  data  sets)  , 

(value  of  the  coupling  coefficient)  , 

(index  of  the  first  mapping  data  set)  , 

(index  of  the  second  mapping  data  set)  , 

(number  of  unknowns  per  reduced  unkn.)  , 

(g-element  — >  unkn.  mapping  index)  , 

(reduced  unknown  index)  , 

(number  of  tensor  multipl.  table)  , 

(first  tensor  index)  , 

(second  tensor  index)  , 

(tensor  index  of  the  basic  matrix  element)  , 
(unknown  number)  . 


=  kkxcj  *  nrj  +  nij 

(7.2) 

Running  summation  (loop)  indices  are  marked  with  boldface  descriptions. 

The  main  point  of  using  the  scheme  of  computation  of  Eq.(7.1)  and  the  code  (7.18)  is  to 
minimize  the  amount  of  work  and  storage  in  the  most  expensive  operations  of  computing 
the  “basic”  matrix  elements  (between  the  “elementary”  basis  functions).  A  set  of  such 
elements,  with  all  required  tensor  indices  I,  is  evaluated  for  the  given  pair  of  g-elements,  q1 
and  72 ,  as  soon  as  these  indices  are  set,  i.e.,  outside  the  sums /loops  over  the  indices  u,  k. 
r  j ,  and  r2.  This  set  of  values  is  only  stored  temporarily,  and  its  contributions  to  all  blocks 
and  elements  of  the  matrix  are  evaluated  in  the  inner  loops  (through  u,  k,  rq,  r2). 
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The  procedure  of  generating  input  data  for  the  computation  of  the  stiffness  matrix 
blocks,  and  the  actual  computation  is  shown,  schematically,  in  Fig.  3.  We  will  repeatedly 
refer  to  this  Figure  and  explain  its  details  in  the  following. 


StorelCOO 

tensor 

multiplication 
tables 


pcsSetGC_... 

direct  mapping  inverse  mapping 

*  unkns.  to  g-elements  g-elements  to  unkns. 


pcot_ 


t 

AddCb 


coupling  data 


pcaSetCl_...  pcaSetC2_ 

mapping  data  l|  mapping  data  2  | 


(ppcb 


PPcal)  ppca2) 


SetAAnCC 

I 

~  matrix  blocks 

(ppcsa 


Figure  3:  Schematics  of  construction  of  matrix  blocks.  Data  are  marked  with  ovals,  the 
remaining  entries  are  routines  called  in  the  code. 

7.2  Data  structures 

Geometry  elements,  basis  functions,  and  unknowns.  As  an  input  to  construction 
of  the  stiffness  matrix,  we  consider  two  types  of  geometries: 

1.  Two-dimensional  manifolds  discretized  with  triangles.  The  triangles  (or  facets,  de¬ 
noted  by  /)  are  referred  to  as  s-elements  (surface  elements)  and  the  edges  as  the  as 
b-elements  (boundary  elements).  Those  geometries  form  boundaries  of  and  inter¬ 
faces  separating  various  material  regions. 

2.  Three-dimensional  manifolds  discretized  with  tetrahedra  (a  generalization  to  hexahe- 
dra  is  possible  but  only  partly  implemented).  In  this  case  the  tetrahedra  (denoted 
by  t )  are  referred  to  as  s-elements  and  the  their  boundaries  (i.e.,  triangular  facets, 
denoted  by  /)  as  b-elements. 

Hence, 


2- dim:  s-element:  facet  /  b-element:  edge  e  , 

3- dim:  s-element:  tetrahedron  t  b-element:  facet  /  . 


(7.3a) 

(7.3b) 
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We  refer  to  all  these  geometry  elements  as  “g-elements” .  The  elements  are  oriented,  i.e. ,  an 
edge  are  defined  as  an  ordered  pair  of  vertices,  and  a  facet  as  an  ordered  triplet,  defining 
the  direction  of  the  normal  according  to  the  right-handed  screw  rule. 

In  the  code  s-elements  are  usually  indexed  by  ns  and  b-elements  by  nb.  The  index  ng 
is  normally  used  for  general  g-elements. 

Geometries  are  described  in  the  code  by  the  structure  CGEO,  whose  main  elements  are: 

1.  cdim:  number  of  dimensions  of  the  geometry  (as  a  manifold),  i.e.,  1  for  a  curve,  2  for 
a  surface,  3  for  a  volume. 

2.  nnv:  number  of  vertices. 

3.  nns:  number  of  s-elements. 

4.  nnb:  number  of  b-elements. 

5.  nsvx:  number  of  vertices  per  s-element,  e.g.,  4  for  tetrahedra. 

6.  nbvx:  number  of  vertices  per  b-element,  e.g.,  3  for  triangles  (facets). 

7.  nnsd:  number  of  sides  of  an  s-element,  e.g.,  4  for  tetrahedra. 

8.  pvx:  array  of  sequentially  stored  (x,  y,  z)  coordinates  of  the  vertices,  hence  3  *  nnv 
floating  point  numbers. 

9.  psv:  array  of  numbers  of  vertices  defining  s-elements,  hence  nsvx  *  nns  integers, 
with  values  from  1. 

10.  pbv:  array  of  numbers  of  vertices  defining  b-elements,  hence  nbvx  *  nnb  integers, 
with  values  from  1. 

11.  psb:  array  of  numbers  of  b-elements  forming  boundaries  of  s-elements,  hence  nnsd  * 
nns  integers,  with  signed  values  from  ±1.  The  signs  define  orientation  of  a  b-element 
with  respect  to  the  s-element.  E.g.,  for  cdim  =  3,  psb  >  0  if  the  facet  (b-element) 
orientation  is  away  from  the  tetrahedron  (s-element).  Similarly,  for  cdim  =  2,  psb  >  0 
if  the  edge  (b-element)  orientation  relative  to  the  facet  (s-element)  is  counter-clockwise 
(looking  from  the  positive  side  of  the  facet). 

12.  pbs:  array  of  numbers  of  s-elements  adjacent  to  b-elements,  hence  2  *  nnb  integers, 
with  values  from  1.  The  first  entry  in  each  pair  is  the  number  of  the  s-element  on  the 
negative  side  of  the  b-element.  If  the  b-element  is  a  facet,  its  negative  side  is  opposite 
the  direction  of  the  normal.  If  it  is  an  edge,  its  negative  side  is  the  left  side  when 
looking  from  the  positive  side  of  the  surface.  If  the  b-element  has  only  one  adjacent 
s-element,  the  missing  entry  in  the  pair  is  zero. 

We  also  consider  several  types  of  “fundamental”  basis  functions  associated  with  s- 
elements  of  the  geometry: 

(i)  For  cdim  =  2: 
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1.  Constant  scalar  basis  functions  supported  on  facets. 

2.  Linear  scalar  basis  functions  supported  on  facets. 

3.  Linear  vector  basis  functions  supported  on  facets. 

(ii)  For  cdim  =  3: 

1.  Constant  scalar  basis  functions  supported  on  tetrahedra. 

2.  Linear  scalar  basis  functions  supported  on  tetrahedra. 

3.  Linear  vector  basis  functions  supported  on  tetrahedra. 

The  above  basis  functions  can  be  combined  to  form  “composite”  basis  functions  sup¬ 
ported  on  sets  of  s-elements,  and  associated  with  various  g-elements.  There  are  also  relations 
between  some  of  these  basis  functions  and  derivatives  of  other  functions. 

Material  parameters.  We  assume  material  parameters  constant  over  s-elements.  In  the 
code  they  are  described  by  the  pair  (kkms,  pms),  where  int  kkms  is  the  number  of  parameters 
per  s-element,  and  COMPLEX  pms  []  is  the  array  of  parameter  values  stored  sequentially  for 
all  s-elements  (i.e. ,  kkms  *  nns  complex  numbers).  For  example,  in  acoustics  we  have 
kkms  =  2,  corresponding,  e.g.,  to  relative  density  and  relative  compressibility  of  the  material. 
In  elasticity  we  need  kkms  =  3  in  order  to  access  relative  density,  and  relative  values  of  Lame 
coefficients. 

Tensor  tables.  The  tables  used  to  combine  tensor  indices  of  basis  functions  into  in¬ 
dices  of  matrix  elements  are  stored  in  predefined  integer  matrices  in  the  coordinate  format, 
MI_C00  *pcot.  They  appear  in  Fig.  3  under  the  generic  name  pcot_.  ...  A  set  of  various 
tables  is  defined  and  then  used  as  input  to  further  operations. 

Unknown  <-»  g-element  mapping  tables.  These  tables  are  stored  as  integer  matrices 
in  the  sparse-row  format,  MI_CSR  *pcsg  (for  the  “forward”  mapping  from  unknowns  to 
g-elements)  and  MI_CSR  *pcsr  (for  the  inverse  mapping).  More  details  on  these  mappings 
are  given  below  in  the  description  of  the  structure  CCFA  of  the  mapping  coefficients. 

Unknown  g-element  mapping  data.  Sets  of  tables  coefficients  relating  basis  func¬ 
tions  to  g-elements  (such  as  ppcal  and  ppca2  in  Fig.  3)  are  stored  in  structures  CCFA  defined 
in  the  header  vmax_ccfa.h  as 
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/*  structure  for  storing  C  coefficients  —  prefix  ca  */ 
typedef  struct 
{ 

char  name [80] ;  /*  name  to  identify  the  object  */ 

int  kxt  [3] ;  /*  unknown  type: 

kxt  [0]  —  unknowns  associated  with: 

O(vertex)  l(edge)  2(facet)  3 (tetrahedron) 
kxt[l]  —  physical  quantity  (tensor  rank) 

O(pressure)  1 (displacement)  3(stress  tensor) 
kxt  [2]  —  physical  problem  type 

*/ 

int  kgt ;  /*  g-element  type: 

O(vertex)  l(edge)  2(facet)  3 (tetrahedron)  */ 
int  kbt ;  /*  basis  function  type:  0 (const)  1 (linear)  */ 

int  nnx;  / *  number  of  unknowns  */ 

int  nnr;  /*  reduced  number  of  unknowns  */ 

int  kkx;  / *  number  of  unknowns  per  "reduced  unknown"  */ 

int  nni;  /*  additional  number  of  tensor  components  (not  used)  */ 

/*  mapping:  unknowns  — >  geometry  elements 

—  allocation:  nnr  (if  nnr  >  0)  or  nnx  */ 

MI_CSR  *pcsg;  /*  arrays  of  lists  of  g-elements  for  (red.)  unkns . 

*  row  #  (from  0)  =  nr  =  reduced  unknown  # 

*  col  #  (from  1)  =  ng  =  g-element  # 

*  value  (from  1)  =  kb  =  #  of  the  basis  fn.  (<  kkx) 

*/ 

/*  inverse  mapping:  geometry  elements  — >  unknowns 

—  allocation:  nng 

constructed  as  transpose  of  pcsg,  except  that  the  values 
are  addresses  nag  in  pcsg  */ 

MI_CSR  *pcsr;  /*  arrays  of  lists  of  (red.)  unknowns  for  g-elems. 

*  row  #  (from  0)  =  ng  =  g-element  # 

*  col  #  (from  1)  =  nr  =  reduced  unknown  # 

*  value  (from  1)  =  nag  =  storage  #  in  pcsg 

*/ 

/*  coefficient  values:  */ 

COMPLEX  *Pgc;  / *  array  of  values  with  access  given  by  pcsg: 

C  =  pgc  [nag]  [kx]  [ni] 

for  each  element  of  pcsg:  [kkx] [nni]  values 
(nni  not  used)  */ 

>  CCFA; 

(7.4) 

In  this  structure: 

1.  The  array  is  given  a  name  name  for  identification. 

2.  The  array  of  three  integers  kxt  [3]  characterize  the  unknown  type:  the  g-element  with 
which  the  unknown  is  associated  and  the  tensor  rank  related  to  the  given  physical 
quantity;  the  third  number  is  currently  not  used  (it  is  set  to  an  undefined  value,  -1). 
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3.  The  integer  kgt  defines  the  g-element  type. 

4.  The  integer  kgt  defines  the  basis  function  type. 

5.  The  integer  nnx  denotes  the  number  of  unknowns. 

6.  The  integer  nnr  denotes  the  “reduced  number  of  unknowns”,  i.e. ,  generally,  the  num¬ 
ber  of  unknowns  not  counting  the  number  of  their  tensor  components.  E.g.,  for  vector 
displacement  unknowns  nnr  is  equal  to  the  number  of  vertices,  while  nnx  is  three 
times  that  number.  By  using  nnr  instead  of  nnx  we  can  avoid  replicating  identical 
coefficients  for  several  tensor  components. 

7.  The  related  entry,  kkx,  is  the  number  of  unknowns  per  reduced  unknown,  e.g.,  3  for 
a  vector  displacement  unknown. 

8.  The  entry  nni  is  the  number  of  tensor  components  characterizing  the  basis  function. 

9.  The  integer  sparse-row  format  matrix  MI_CSR  *pcsg  defines  the  mapping  from  re¬ 
duced  unknowns  to  geometry  elements:  its  rows  are  indexed  by  reduced  unknown 
numbers,  columns  by  g-element  numbers,  and  the  values  store  the  numbers  of  basis 
functions  (up  to  kkx).  Usually,  to  one  unknown  there  correspond  several  g-elements; 
e.g.,  a  vector  displacement  unknown  associated  with  a  vertex  is  supported  by  a  set  of 
tetrahedra  sharing  that  vertex. 

10.  The  integer  sparse-row  format  matrix  MI_CSR  *pcsr  defines  the  mapping  inverse  to 
the  previous  one,  i.e.,  from  geometry  elements  to  the  reduced  unknowns:  its  rows 
are  indexed  by  g-element  numbers,  and  columns  by  reduced  unknown  numbers.  In 
this  case  the  values  stored  in  the  the  matrix  are  storage  numbers  of  elements  in  the 
previous  matrix,  pcsg.  This  information  is  used  in  constructing  matrix  elements  in  a 
loop  through  g-elements  in  the  routine  SetAAnCC,  as  described  below. 

11.  Finally,  the  coefficient  values  are  stored  in  the  complex  array  pgc.  The  array  contains 
kkx  *  nni  values  for  each  element  of  the  mapping  array  pcsg.  These  values  may,  in 
general,  depend  on  the  material  parameters. 

Tables  of  coupling  coefficients.  Each  MoM  matrix  block6  is  constructed  by  specifying 
a  set  of  “coupling  coefficients”  which  specify  how  to  combine  products  of  basis  functions 
(associated  with  unknowns)  into  matrix  elements.  These  coefficients  are  stored  in  a  structure 

6Compressed  matrix  blocks  will  be  created  in  a  similar  scheme. 
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/*  structure  for  storing  coefficients  B  and  multiplication  tables  T 
—  prefix  cb  */ 
typedef  struct 
{ 

char  name  [80]  ;  / *  name  to  identify  the  object  (matrix  block)  */ 

MC_C00  *pcob;  / *  array  of  coefficients  B  defining  couplings 

*  row  #  (from  0)  =  nul  =  number  of  the  set  C~1 

*  col  #  (from  1)  =  nu2  =  number  of  the  set  C~2 

*  value  =  b  =  complex  value  of  the  coef f . 

*  elements  are  assumed  to  be  sorted  by  rows  and  cols 

*/ 

MI_C00  **ppcot;  /*  array  of  index  (tensor)  multiplication  tables 

*  one  MI_C00  matrix  for  each  element  of  pcob 
*/ 

TNAME  *ptname ;  / *  array  of  names  of  multiplication  tables 

*  one  name  for  each  element  of  pcob 
*/ 

>  CBFA; 

(7.5) 

In  this  structure: 

1.  The  name  name  identifies  the  set  of  coupling  coefficients  and  thus  the  resulting  matrix 
block. 

2.  The  complex  coordinate-format  matrix  MC_C00  *pcob  specifies  a  set  of  pairs  of  un¬ 
known  4-4-  g-element  mapping  data  (C^1),  C^),  identified  by  their  numbers,  and  com¬ 
plex  coefficients  with  which  contributions  of  those  pairs  are  to  be  added  to  the  result. 

3.  To  each  element  of  the  matrix  MC_C00  *pcob  there  corresponds  a  pointer  to  a  tensor 
multiplication  table. 

4.  Names  of  the  above  multiplication  tables  are  also  stored  in  an  array. 

7.3  Operations 

We  describe  here  the  main  operations  executed  by  the  code  in  assembling  blocks  of  MoM 
matrix  elements,  according  to  the  procedure  represented  in  Fig.  3. 

Structure  of  the  matrix  construction  “driver”.  We  reproduce  below  parts  of  the 
present  main  program  which  serves  as  a  driver  for  generating  blocks  of  the  MoM  matrix. 
This  code  is  a  simple  illustration  of  the  data  flow  shown  in  Fig.  3. 
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//  set  up  material  data:  kkmsl,  pmsl,  kkms2,  pms2 

//  set  up  two  tensor  multiplication  tables 
MI_C00  *pcot  =  NULL; 

//  s  s  ->  s  (scalar  scalar  ->  scalar) 

StorelCOO (pcot_sss ,  1,  1,  1); 

//  v  v  ->  t  (vector  vector  ->  symmetric_tensor) 

StorelCOO (pcot_vvt ,  1,  1,  1);  StorelCOO (pcot_vvt ,  1,  2,  2); 

StorelCOO (pcot_vvt ,  1,  3,  3); 

StorelCOO (pcot_vvt ,  2,  2,  4);  StorelCOO (pcot_vvt ,  2,  3,  5); 

StorelCOO (pcot_vvt ,  3,  3,  6); 

//  create  nnu  =  2  sets  of  coupling  data  ppcb[] 
nnu  =  2 ; 

//  first  set 
pcb  =  ppcb  [0] ; 

AddCb(pcb,  1,  1,  C0MPLEX(1.),  pcot_sss,  "t_sss"); 

AddCb(pcb,  2,  1,  C0MPLEX(2.),  pcot_sss,  "t_sss"); 

//  second  set 
pcb  =  ppcb  [1] ; 

AddCb(pcb,  1,  1,  C0MPLEX(2.),  pcot_vtv,  "t_vtv"); 

AddCb(pcb,  2,  1,  C0MPLEX(1.),  pcot_vtv,  "t_vtv"); 

//  set  up  3  unknown  -  g-element  mapping  tables  pcsg. . . 

//  and  inverse  mapping  tables  pcsr. . . 

//  for  Gl: 

pcsgC_f_dTvD  =  pcsSetGC_f _dTvD(pcgeol) ; 
pcsrC_f_dTvD  =  pcsSetR(pcsgC_f _dTvD) ; 
pcsgL_f_Fv  =  pcsSetGL_f _Fv(pcgeol) ; 
pcsrL_f _Fv  =  pcsSetR(pcsgL_f _Fv) ; 

//  for  G2 : 

pcsgC_t_Tv  =  pcsSetGC_t_Tv(pcgeo2) ; 
pcsrC_t_Tv  =  pcsSetR(pcsgC_t_Tv) ; 

//  create  nncal  =  2  sets  of  mapping  data  for  Gl,  nnca2  =  1  sets  for  G2 
nncal  =  2;  nnca2  =  1; 

//  create  elements  of  arrays  of  mapping  data 
//  for  Gl: 

ppcal [0]  =  pcaSetCl_C_f (pcgeol ,  kkmsl,  pmsl, 

pcsgC_f _dTvD,  pcsrC_f _dTvD ,  "C_f"); 
ppcal [1]  =  pcaSetCl_L_f (pcgeol ,  kkmsl,  pmsl, 

pcsgL_f_Fv,  pcsrL_f_Fv,  "L_f"); 

//  for  G2 : 

ppca2 [0]  =  pcaSetC2_C_t (pcgeo2 ,  kkms2,  pms2, 

pcsgC_t_Tv,  pcsrC_t_Tv,  "C_t"); 

//  compute  nnu  =  2  matrix  blocks  ppcsAn[] 

SetAAnCC(akO,  pcgeol,  pcgeo2,  nnu,  ppcb,  nncal,  nnca2,  ppcal,  ppca2,  &ppcsAn) ; 

(7.6) 
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Construction  of  tensor  multiplication  tables.  In  the  initial  stage  of  code  execution 
we  define  a  set  of  tables  for  operations  on  tensor  indices  of  the  basis  functions  and  matrix 
elements.  They  are  generated  “by  hand”  by  calling  the  routine  StorelCOO  (Fig.  3)  and 
specifying  nonzero  elements  of  the  tables.  E.g.,  the  table  specifying  a  symmetric  tensor  in 
terms  of  two  vectors  is  constructed  by 

//  v  v  ->  t  (vector  vector  ->  symmetric_tensor) 
pcot  =  pcoICOOzO; 

StorelCOO (pcot ,  1,  1,  1); 

StorelCOO (pcot ,  2,  2,  4); 

StorelCOO (pcot ,  2,  3,  5); 

StorelCOO (pcot ,  1,  2,  2); 

StorelCOO (pcot ,  1,  3,  3); 

StorelCOO (pcot ,  3,  3,  6); 

MI_C00  *pcot_vvt  =  pcoSortICOO (pcot) ;  (7.7) 

These  instructions  (in  which  table  elements  may  be  entered  in  any  order  -  they  are  sorted 
in  the  result)  create  a  table  named  pcot_vvt  which  is  shown  in  the  test  output  as 

T_vvt : 

12  3 


112  3 

2  4  5 

3  6  (7.8) 

It  indicates  how  the  pairs  of  components  of  two  vectors  are  assigned  to  the  six  independent 
components  of  the  symmetric  tensor;  the  row  numbers  refer  to  the  first,  and  the  column 
numbers  to  the  second  vector. 

Construction  of  mapping  tables.  The  code  constructs  a  set  of  tables  describing  map¬ 
ping  between  the  unknowns  and  the  g-elements  (pcsgC_.  .  .  and  pcsgC_.  .  .  in  Fig.  3). 

The  forward  mappings  are  generated,  for  specific  cases,  by  separate  routines,  called  with 
one  of  the  geometries  as  the  argument.  For  example,  the  forward  mapping  for  linear  basis 
functions  supported  on  faces  /  in  sets  of  faces  J~v  sharing  a  vertex  v,  in  the  first  geometry, 
is  generated  by  the  routine 

//  Returns  pcsg  for  linear  b.f.,  f  in  F_v. 

/ /  input : 

//  pcgeo  =  geometry 

MI_CSR  *pcsSetGL_f _Fv (CGEO  *pcgeo)  (7.9) 

in  vmax_cf  f  a_c  .  cpp.  It  is  called  as 

MI_CSR  *pcsgL_f_Fv  =  pcsSetGL_f _Fv(pcgeol) ;  (7.10) 
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and  stores  the  output  as  the  array  pcsgL_f  _Fv. 

All  inverse  mappings  are  generated  from  the  forward  ones  by  calling  the  same  routine 
pcsSetR.  E.g.,  in  the  considered  example,  the  inverse  mapping  is  obtained  by 

pcsrL_f_Fv  =  pcsSetR(pcsgL_f _Fv) ;  (7-11) 


Construction  of  mapping  data.  The  final  set  of  mapping  data,  representing  relations 
between  the  unknowns  and  g-elements  (ppcal  and  ppca2  in  Fig.  3) 

Specific  tables  of  coefficients  C  for  various  sets  of  unknowns  and  basis  functions  are 
generated  by  specialized  routines.  For  example,  a  set  of  coefficients  C W  for  linear  basis 
functions  supported  on  faces  is  generated  by  the  routine 


//  Computes  and  returns  a  set  of  C  coefficients  for  a  specific  set 
//  of  unknowns  and  basis  functions 
/ /  input : 

=  geometry 

=  number  of  material  parameters  per  s-element 
=  array  of  material  parameters  values  for  s-elements 
=  a  mapping  table:  unknowns  — >  g-elements 
=  a  mapping  table:  g-elements  — >  unknowns 
=  name  of  the  set 


//  pcgeo 
//  kkms 
//  pms 
//  pcsg 
//  pcsr 
//  szName 
// 


//  C  for  set  1  of  linear  basis  functions  supported  on  facets 
CCFA  *pcaSetCl_C_f (CGEO  *pcgeo,  int  kkms,  COMPLEX  *pms, 

MI_CSR  *pcsg,  MI_CSR  *pcsr,  char  *szName) 


(7.12) 


It  may  be  called  with  various  unknown  g-element  mappings.  E.g.,  the  call 


peal  =  pcaSetCl_L_f (pcgeol ,  kkmsl,  pmsl, 

pcsgL_f_Fv,  pcsrL_f_Fv,  "L_f");  (7-13) 

creates  the  coefficient  set  named  L_f  by  using  as  arguments  the  geometry  pcgeol,  material 
parameters  specified  by  kkmsl  and  pmsl,  and  the  mapping  tables  pcsgL_f  _Fv  and  its  inverse 
pcsrL_f_Fv.  These  mapping  tables  have  been  created  previously  (Eqs.  (7.10)  and  (7.11)) 
for  basis  functions  supported  on  faces  /  in  sets  of  faces  J-v  sharing  a  vertex  v. 

One  of  the  sets  of  the  mapping  data  is  shown  in  the  test  output  as 
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C~ (1)2: 


name :  L_f 


kxt :  0  1  -1 

kgt:  : 

2 

kbt : 

1 

nnx : 

375 

nnr : 

125 

kkx : 

3 

nni : 

3 

pcsg: 

125  rows, 

1020  columns,  3060  elements 

rows :  1  . 

.  .  125  (125),  columns:  1  ...  1020  (1020) 

pcsr : 

1020  rows 

,  125  columns,  3060  elements 

rows :  1  . 

.  .  1020  (1020),  columns:  1  ...  125  (125) 

#g: 

3060  #g  /  #r:  24.5 

pgc :  27540  elements  |C|:  13139.8  (7.14) 


Construction  of  coupling  data.  The  coupling  data  (an  array  of  which  is  denoted  by 
ppcb  in  Fig.  3)  are  constructed  “by  hand”  by  selecting  specific  sets  of  mapping  data  pairs 
(C^1),  C^),  the  corresponding  coefficients,  and  the  tensor  multiplication  tables.  For  in¬ 
stance,  the  second  set  of  coupling  data  in  the  example  code  is  generated  by  the  instructions 

pcb  =  ppcb  [1] ; 

AddCb(pcb,  1,  1,  C0MPLEX(2.),  pcot_vtv,  "t_vtv"); 

AddCb(pcb,  2,  1,  COMPLEX (1.),  pcot_vtv,  "t_vtv") ;  (7.15) 

The  result  of  these  assignments  is  shown  in  the  test  output  as 

B~2_mil_nu2 


nul 

nu2 

Re  B 

Im  B 

C~(l) 

C~(2) 

T 

1 

1 

2 

0 

O 

1 

Hh 

C_t 

t_vtv 

2 

1 

1 

0 

L_f 

C_t 

t_vtv 

(7.16) 


Construction  of  a  set  of  matrix  blocks.  The  final  stage  of  assembling  a  set  of  nnu 
matrix  blocks  ppcsa  (Fig.  3)  is  carried  out  by  the  routine 
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//  Allocates  and  computes  set  of  near-field  matrices  ppcsa 
//  for  a  set  of  coupling  coefficients  B  =  ppcb, 

//  and  mapping  coefficients  ppcal  and  ppca2 
//  input : 


//  akO 
//  pcgeol 
//  pcgeo2 
//  nnu 
//  ppcb 
//  nncal 
//  nnca2 
//  ppcal 
//  ppca2 
//  output : 

//  ppcsa  =  array  of  nnu  matrix  blocks 

void  SetAAnCC (FL0AT_T  akO,  CGEO  *pcgeol,  CGEO  *pcgeo2, 
int  nnu,  CBFA  **ppcb, 

int  nncal,  int  nnca2,  CCFA  **ppcal,  CCFA  **ppca2, 
MC_CSR  ***pppcsa) 


wave  number 
geometry  1 
geometry  2 

number  of  matrices  to  be  created 
array  of  nnu  coupling  coefficients  B 
number  of  sets  of  mapping  coefficients 
number  of  sets  of  mapping  coefficients 
array  of  nncal  sets  of  coefficients  C~1 
array  of  nnca2  sets  of  coefficients  C~2 


(7.17) 


in  vmax_ccf a. cpp. 

The  routine  uses  as  input: 

1.  The  wave  number  akO. 

2.  The  geometries  G1  and  G2  (pcgeol  and  pcgeo2),  set  of  nnu  coupling-coefficient  data 

ppcb. 

3.  Numbers  nncal  and  nnca2  of  sets  of  mapping  coefficients  C ^  and  C ^  to  be  used  in 
the  construction  (they  are  only  used  to  check  index  values). 

4.  Arrays  ppcal  and  ppca2  of  mapping  coefficient  data  C ^  and  C^2\ 

The  output  generated  by  the  routine  -  a  set  of  nnu  matrices  -  is  stored  in  the  array  of 
matrices  ppcsa. 

The  general  structure  of  the  routine  is  summarized  in  the  pseudo-code 
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//  loops  through  g-element  types  for  geometries  G1  and  G2 
for  kgtl  =  1,  ...  ,3 

nngl  =  nngklfkgtl]  //  number  of  g-elements  of  type  kgtl 
//  get  pointer  to  the  required  type  kgtl  of  g-elements  for  geom.  Gl: 
pgl  =  pcgeol->  ■  ■  ■ 
for  kgt2  =  1,  ...  ,3 

nngl  =  nngklfkgtl]  //  number  of  g-elements  of  type  kgtl 

//  get  pointer  to  the  required  type  kgt2  of  g-elements  for  geom.  G2: 

Pg2  =  pcgeo2->  ■  ■  ■ 

//  loops  through  g-elements  in  geometries  Gl  and  G2 
for  ngl  =  1,  . .  .  ,  nngl 

for  ng2  =  1,  ...  ,nng2 

/ /  compute  set  of  basic  matrix  elements,  store  in  paan 
SetAglg2(pgl,  ngl,  pg2,  ng2,  paan) 

/ /  loop  through  numbers  nu  of  blocks  to  be  generated 
for  nu  =  1,  ...  ,  nnu 

peb  =  ppcb  [nu]  //  coupling  coefficient  data  for  the  matrix  block  nu 

pcob  =  pcb->pcob  //  array  of  coupling  coefficients  for  sets  of  basis  functions 

kkee  =  pcob->na  //  number  of  pairs  of  mapping  data  sets 

ppcot  =  pcb->ppcot  //  set  of  tensor-component  multiplication  tables 

//  loop  through  pairs  of  mapping  data  sets  ( C('1\C<'2'1 ) 

for  kcc  =  1,  . . .  ,  kkee 

//  get  coupling  coefficient  B  and  tensor  multiplication  table 
//  (for  later  use) 

cb  =  pcob->paa[kcc]  //  coupling  value 

pcot  =  ppcot  [kcc]  //  tensor  multiplication  table 

nnat  =  nalCOO(pcot)  //  number  of  its  elements 

//  get  numbers  of  mapping  data  sets  (C^,  C^2>)  (row  and  col.) 

ncal  =  pcob->pia[kcc] 

nca2  =  pcob->pja[kcc] 

peal  =  ppcal  [ncal]  //  mapping  data  1 

pca2  =  ppca2  [nca2]  //  mapping  data  2 

/  /  get  numbers  of  unknowns  per  reduced  unknown 

kkxcl  =  pcal->kkx 

kkxc2  =  pca2->kkx 

pcsrl  =  pcal->pcsr  //  inverse  mapping  table  1 

pcsr2  =  pca2->pcsr  //  inverse  mapping  table  2 

pgcl  =  pcal->pgc  //  mapping  coefficients  1 

pgc2  =  pca2->pgc  //  mapping  coefficients  2 

//  get  bounds  for  loops  through  inverse  mappings  elements 

nnarll  =  pcsrl->pia[ngl] 

nnarl2  =  pcsrl->pia[ngl  +  1] 

nnar21  =  pcsr2->pia[ng2] 

nnar22  =  pcsr2->pia[ng2  +  1] 
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/ /  loop  through  elements  of  of  the  inverse  map  1 
for  narl  =  nnarll,  .  .  .  ,nnarl2 

nrl  =  pcsrl->pja[narl]  //  reduced  unkn.  1  number 
nagl  =  pcsrl->paa[narl]  //  dir.  map.  storage  number  1 
/ /  loop  through  elements  of  of  the  inverse  map  2 
for  nar2  =  nnar21,  ....  nnar22 

nr2  =  pcsr2->pja[nar2]  //  reduced  unkn.  2  number 
nag2  =  pcsr2->paa[nar2]  //  dir.  map.  storage  number  2 
//  loop  through  elements  of  the  tensor  multipl.  table 
for  nat  =  1 ,  ...  ,  nnat 

//  get  tensor  indices 
nil  =  pcot->pia[nat] 
ni2  =  pcot->pja[nat] 

//  get  basic  matrix  element  tensor  index 
nia  =  pcot->paa[nat] 

//  get  unknown  numbers 
nxl  =  kkxcl  *  nrl  +  nil 

(7.18) 

nx2  =  kkxc2  *  nr2  +  ni2 
//  get  mapping  coefficients 
acl  =  pgcl [kkxcl  *  nagl  +  nil] 
ac2  =  pgc2 [kkxc2  *  nag2  +  ni2] 

//  multiply  basic  matrix  element  by  coefficients 
cba  =  cb  *  acl  *  ac2  *  paanfnia] 

//  add  to  element  in  output  block  nu 
StoreAddCSRna(ppcsa[nu] ,  nxl,  nx2,  cba) 
endfor  nat 
endfor  nar2 
endfor  narl 
endfor  kcc 
endfor  nu 
endfor  ng2 
endfor  ngl 
endfor  kgt2 
endfor  kgtl 


8  Summary  of  the  developments  in  formulation  and  imple¬ 
mentation  of  integral  equations  for  elasticity 

We  have  described  the  formulation  of  several  types  of  integral  equations  we  are  implementing 
in  our  solver,  and  the  detailed  expressions  for  the  resulting  matrix  equations.  As  the 
presented  material  shows,  the  structure  of  the  equations  is  quite  complex  and  requires  a 
significant  amount  of  analytic  and  programming  work. 
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In  the  area  of  code  implementation,  we  have  designed  a  general  and  flexible  solution 
scheme  involving  construction  of  the  stiffness  matrix,  its  compression,  and  matrix-vector 
multiplication  routines,  which  allows  relatively  easy  implementation  of  various  types  of 
integral  equations,  including  volumetric  and  surface  equations,  equations  based  on  first-  and 
second-order  formulations,  and  equations  specially  adapted  to  high-contrast  problems,  hence 
involving  different  unknown  fields  and  treating  material  properties  in  various  ways.  Our 
general  goal,  however,  is  not  to  unnecessarily  restrict  the  allowable  types  of  equations,  but 
rather  keep  the  scheme  open-ended  and  allow  incorporating  new  formulations  in  a  possibly 
straightforward  way.  To  this  end,  we  are  developing  a  comprehensible  and  extensible  library 
of  routines  for  constructing  particular  blocks  and  sets  of  matrix  elements,  as  well  as  their 
compressed  representations,  appearing  in  various  integral-equation  formulations.  These 
routines  create  input  data  used  by  a  general  routine  whose  task  is  to  assemble  the  entire 
matrix  and  store  it  in  a  compressed  form;  the  overall  scheme  of  the  matrix  construction  is 
visualized  in  Fig.  3  and  exemplified  with  a  set  of  routines  described  in  Sec.  7.3. 
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